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Abstract 

The concept of the limiting step gives the hmit simplification: the whole network behaves as a single step. This is 
the most popular approach for model simplification in chemical kinetics. However, in its elementary form this idea is 
applicable only to the simplest linear cycles in steady states. For such simple cycles the nonstationary behaviour is 
also limited by a single step, but not the same step that limits the stationary rate. In this paper, we develop a general 
theory of static and dynamic limitation for all linear multiscale networks. Our main mathematical tools are auxiliary 
discrete dynamical systems on finite sets and specially developed algorithms of "cycles surgery" for reaction graphs. 
New estimates of eigenvectors for diagonally dominant matrices are used. 

Multiscale ensembles of reaction networks with well separated constants are introduced and typical properties of 
such systems are studied. For any given ordering of reaction rate constants the explicit approximation of steady 
state, relaxation spectrum and related eigenvectors ("modes") is presented. In particular, we prove that for systems 
with well separated constants eigenvalues are real (damped oscillations are improbable). For systems with modular 
structure, we propose the selection of such modules that it is possible to solve the kinetic equation for every module 
in the explicit form. All such "solvable" networks are described. The obtained multiscale approximations, that we 
call "dominant systems" are computationally cheap and robust. These dominant systems can be used for direct 
computation of steady states and relaxation dynamics, especially when kinetic information is incomplete, for design 
of experiments and mining of experimental data, and could serve as a robust first approximation in perturbation 
theory or for preconditioning. 

Key words: Reaction network, multiscale ensemble, dominant system, multiscale asymptotic, model reduction, spectral gap 
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1. Introduction populi, vox Dei," especially in the epoch of citation 

indexes, impact factors and bibliometrics. Let us ask 
Which approach to model reduction is the most Google. It gave on 31st December 2006: 
important? Population is not the uhimate judge, " for "quasi-equilibrium" - 301000 links; 
and popularity is not a scientific criterion, but "Vox " for "quasi steady state" 347000 and for "pseudo 



steady state" 76200, 423000 to gether; 

for ou r f avorite "slow manif ol d" (TCorban fc Karlin 
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98100; 

- for such a framework topic as "singular perturba- 
tion" Google gave 361000 links; 

- for "model reduction" even more, as we did ex- 
pect, 373000; 

- but for "limiting step" almost two times more - 
714000! 

Our goal is the general theory of static and dy- 
namic limitation for multiscale networks. The con- 
cept of the limiting step gives, in some sense, the 
limit simplification: the whole network behaves as a 
single step. As the first result of our paper we intro- 
duce further detail in this idea: the whole network 
behaves as a single step in statics, and as another 
single step in dynamics: even for simplest cycles the 
stationary rate and the relaxation time to this sta- 
tionary rate are limited by different reaction steps, 
and we describe how to find these steps. 

The concept of limitation is very attractive both 
for theorists and experimentalists. It is very useful 
to find conditions when a selected reaction step be- 
comes the limiting step. We can change conditions 
and study the network experimentally, step by step. 
It is very convenient to model a system with lim- 
iting steps: the model is extremely simple and can 
serve as a very elementary building block for further 
study of more complex systems, a typical situation 
both in industry and in systems biology. 

In the lUPAC Compendium of Chemical Termi- 
nology (2007) one can find two articles with a defi- 
ni tion of limitation. 

- Rate-determining sted (2007): Rate-determining 



step (rate-limiting step): "These terms are best 
regarded as synonymous with rate-controlling 

step." 

- iRate-controlling stepi ( 2007 ): "A rate-controlling 
(rate-determining or rate-limiting) step in a reac- 
tion occurring by a composite reaction sequence 
is an elementary reaction the rate constant for 
which exerts a strong effect - stronger than that 
of any other rate constant - on the overall rate." 
It is not wise to object to a definition and here we 
do not object, but, rather, complement the defini- 
tion by additional comments. The main comment is 
that usually when people are talking about limita- 
tion they expect significantly more: there exists a 
rate constant which exerts such a strong effect on 
the overall rate that the effect of all other rate con- 
stants together is significantly smaller. Of course, 
this is not yet a formal definition, and should be 
complemented by a definition of "effect" , for exam- 
ple, by "control function" identified by derivatives 



of the overall rate of reactio n, or by other overall 
rate " sensitivity parameters" ( Rate-controlling sted 
(|2007l )). 

For the lUPAC Compendium definition a rate- 
controlling step always exists, because among the 
control functions generically exists the biggest one. 
On the contrary, for the notion of limitation that is 
used in practice, there exists a difference between 
systems with limitation and systems without limi- 
tation. 

An additional problem arises: are systems with- 
out limitation rare or should they be treated equi- 
tably with limitation cases? The arguments in favor 
of limitation typicality are as follows: the real chem- 
ical networks are very multiscale with very different 
constants and concentrations. For such systems it 
is improbable to meet a situation with compatible 
effects of all different stages. Of course, these argu- 
ments are statistical and apply to generic systems 
from special ensembles. 

During last century, the concept of the limiting 
step was revised several times. First simple idea of 
a "narrow place" (a least conductive step) could be 
applied without adaptation only to a simple cycle 
of irreversible steps t hat are of the fi rst order (see 
C hap. 16 of th e book JohnstoiJ ( 1966h or the paper 
of Bovdl ( 1978f l). When researchers try to apply this 
idea in more general situations they meet various 
difficulties such as: 

- Some reactions have to be "pseudomonomolecu- 
lar." Their constants depend on concentrations of 
outer components, and are constant only under 
condition that these outer components are present 
in constant concentrations, or change sufficiently 
slow. For example, the simplest Michaelis-Menten 
enzymatic reaction is E + S ^ ES E + P {E 
here stands for enzyme, S for substrate, and P 
for product), and the linear catalytic cycle here is 

5 ES S. Hence, in general we must consider 
nonlinear systems. 

- Even under fixed outer components concentra- 
tion, the simple "narrow place" behaviour could 
be spoiled by branching or by reverse reactions. 
For such reaction systems definition of a limiting 
step simply as a step with the smallest constant 
does not work. The simplest example is given by 
the cycle: Ai ^ A2 ^ A3 ^ Ai. Even if the con- 
stant of the last step A3 Ai is the smallest one, 
the stationary rate may be much smaller than k^b 
(where b is the overall balance of concentrations, 

6 = ci -I- C2 -I- C3), if the constant of the reverse 
reaction A2 Ai is sufficiently big. 
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In a series of papers, iNorthrod (|l98ll 120011) 



clearly explained these difficulties with many ex- 
amples based on the isotope effect analysis and 
suggested that the concept of rate-limiting step is 
"outmoded" . Nevertheless, the main idea of limiting 
is so attractive that Northrop's arguments stimu- 
lated the search for modification and improvement 
o f the main c oncept . 

Rav ( 19831 ) proposed the use of sensitivity analy- 



sis. He considered cycles of reversible reactions and 
suggested a definition: The rate-limiting step in a 
reaction sequence is that forward step for which a 
change of its rate constant produces the largest effect 
on the overall rate. In his formal definition of sensi- 
tivity functions the reciprocal reaction rate (l/W) 
and rate constants (l/fcj) were used and the con- 
nection between forward and reverse step constants 

(the equilibrium constant) was ke pt fixed. 

Ray 's approach was revised bv iBrown fc Coopei 
19931 ) from the syst em control analysis point of 



view ( see the book of Cornish-Bowden fc Cardena3 
(|l99nf )'l. They stress again that there is no unique 
rate-limiting step specific for an enzyme, and this 
step, even if it exists, depends on substrate, prod- 
uct and effector concentrations. They demonstrated 
also that the control coefficients 



w 



\W dkj 



[S],[P]. 



where W is the stationary reaction rate and hi are 
constants, are additive and obey the summation 
theorems (as concentrations do). A simple rela- 
tion between control coefficients of rate constants 
and in termediate concentrations was reported by 



Kholodenk o. Westerhoff. fc Brown (1994) . This re- 



lation connects two type of experiments: measure- 
ment of intermediate levels and steady-state rate 
measurements. 

For the analysis of nonlinear cycles the new 



In our approach, we analyze not only the steady 
state reaction rates, but also the relaxation dynam- 
ics of multiscale systems. We focused mostly on the 
case when all the elementary processes have signifi- 
cantly different time scales. In this case, we obtain 
"limit simplification" of the model: all stationary 
states and relaxation processes could be analyzed 
"to the very end" , by straightforward computations, 
mostly analytically. Chemical kinetics is an inex- 
haustible source of examples of multiscale systems 
for analysis. It is not surprising that many ideas 
and methods for such analysis were first invented for 
chemical systems. 

In Sec. [2] we analyze a simple example and the 
source of most generalizations, the catalytic cycle, 
and demonstrate the main notions on this example. 
This analysis is quite elementary, but includes many 
ideas elaborated in full in subsequent sections. 

There exist several estimates for relaxation time 
in chemical reactions (developed, for example, by 
Cheresiz fc Yablonskiil|l983) ). but even for the sim- 



plest cycle with limitation the main property of re- 
laxation time is not widely known. For a simple ir- 
reversible catalytic cycle with limiting step the sta- 
tionary rate is controlled by the smallest constant, 
but the relaxation time is determined by the second 
in order constant. Hence, if in the stationary rate ex- 
periments for that cycle we mostly extract the small- 
est constant, in relaxation experiments another, the 
second in order constant will be observed. 

It is also proven that for cycles with well separated 
constants damped oscillations are impossible, and 
spectrum of the matrix of kinetic coefficients is real. 
For general reaction networks with well separated 
constants this property is proven in Sec.[4l 

Another general effect observed for a cycle is ro- 
bustness of stationary rate and relaxation time. For 
multiscale systems with random constants, the stan- 
dard deviation of constants that determine station- 
ary rate (the smallest constant for a cycle) or re- 
laxation time (the second in order constant) is ap- 



concept of kinetic p olynomial wa s developed 

( Yablonskii. Lazman. fc B vkov (1982i): lLazman fc Yablonpkdj cimatelv n times smaller than the standard de- 
( 199ll )). It was proven that the stationary state of viation of the individual constants (where n is the 



the single-route reaction mechanism of catalytic 
reaction can be described by a single polynomial 
equation for the reaction rate. The roots of the 
kinetic polynomial are the values of the reaction 
rate in the steady state. For a system with limiting 
step the kinetic polynomial can be approximately 
solved and the reaction rate found in the form of 
a series in powers of the l imiting step constant 
(|Lazman fc Yablonskiil (|l988h ). 



cycle length) . Here we deal with the so-called "order 
statistics". This decrease of the deviation as is 
much faster than for the standard error summation, 
where it decreases with increasing n as n~^^^. 

In more general settings, robustness of the re- 
laxati on time was studied by iGorban fc Radulescu 
(|2007| ) for chemical k inetics mod els of genetic and 
signalling networks. iGorban fc R adulescu (2007|) 
proved that for large multiscale systems with hi- 
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erarchical distribution of time scales the variance 
of the inverse relaxation time (as well as the vari- 
ance of the stationary rate) is much lower than the 
variance of the separate constants. Moreover, it can 
tend to faster than 1/n, where n is the number 
of reactions. It was demonstrated that similar phe- 
nomena are valid in the nonlinear case as well. As 
a numerical illustration we used a model of a sig- 
nalling network that can be applied to important 
transcription factors such as NFkB. 

Each multiscale system is characterized by its 
structure (the system of elementary processes) and 
by the rate constants of these processes. To make 
any general statement about such systems when the 
structure is given but the constants are unknown 
it is useful to take the constant set as random and 
independent. But it is not obvious how to chose the 
random distribution. The usual idea to take normal 
or miiform distribution meets obvious difficulties, 
the time scales are not sufficiently well separated. 

The statistical approach to chemical kinetic s 
was developed by iLi. Rosenthal, fc Rabitd (|200lf ): 
Li. Wang. Rabitz. Wang, fc Jaffel (|2002'). and high- 
dimensional model representations (HDMR) were 
proposed as efficient tools to provide a fully 
global statistical analysis of a mode l. The work o f 
Feng. Hooshang i. Chen . Li. Weiss, fc Rabitd (|2004l ) 



was focused on how the network properties are af- 
fected by random rate constant changes. The rate 
constants were transformed to a logarithmic scale 
to ensure an even distribution over the large space. 

The log-uniform distribution on sufficiently wide 
interval helps us to improve the situation, indeed, 
but a couple of extra parameters appears: a = 
minlogfc and /? = max log fc. We have to study 
the asymptotics a — s- — oo, /3 — > oo. This approach 
could be formalized by means of the uniform invari- 
ant distributions of log fc on M". These distributions 
are finite-additive, but not countable-additive (not 
(T-additive). 

The probability and measure theory without 
countable additivity has a long history. In Euclid's 
time only arguments based on finite-additive prop- 
erties of volume were legal. Euclid meant by equal 
area the scissors congruent area. Two polyhedra 
are scissors-congruent if one of them can be cut 
into finitely many polyhedral pieces which can be 
re-assembled to yield the second. But all proofs of 
the formula for the volume of a pyramid involve 
some form of limiting process. Hilbert asked in his 
third problem: are two Euclidean polyhedra of the 
same volume scissors congruent? The answer is 



"n o" (a review of ol d and recent results is presented 
by Neumann ( 19981 )). There is another invariant of 
cutting and gluing polyhedra. 

Finite-additive invariant mea sures on non- 
compact groups wer e studied by Birkhofi ( 1936h 
(see also the book of Hewitt fc Rossi (|l963 ). Chap. 
4). The frequency-based M ises approa c h to prob- 
ability theory foundations (Ivon Mised (Il964l)'l. as 
well a s logical foundations of probabilitv bvl Carnapl 
do not need cr-additivity. Non-Kolmogorov 
probability theories are discussed now in t he con- 
text of quantum ph ysics ( Khrennikov ( 2002f) ). non- 
standard analysis ( Loeb ( 1975f )) and many other 
problems (and we do not pretend provide here a full 
review of related works) . 

We answer the question: What does it mean "to 
pick a multiscale system at random" ? We introduce 
and analyze a notion of multiscale ensemble of reac- 
tion systems. These ensembles with well separated 
variables are presented in Sec. [3l 

The best geometric example that helps us to un- 
derstand this problem is one of Le wis Carroll's P il- 
low Problems pubfished in 1883 (jCarroU l (1958')): 
"Three points are taken at random on an infinite 
plane. Find the chance of their being the vertices of 
an obtuse-angled triangle." (In an acute-angled tri- 
angle all angles are comparable, in an obtuse-angled 
triangle the obtuse angle is bigger than others and 
could be much bigger.) The solution of this prob- 
lem depends significantly on the ensemble definition. 
What does it mean "points are taken at random on 
an infinite plane" ? Our intuition requires translation 
invariance, but the normalized translation invariant 
measure on the plain could not be cr-additive. Nev- 
ertheless, there exist finite-additive invariant mea- 



sures. 



Lewis Carroll proposed a solution that did not 
satisfy some of modern scientists. There exists a 



Guv^( 


1993): Portnovl(ll994l): 


Eisenberg fc Sullivan 


1996h: Falk & Samuel-Cahn 


(200ll)): reduction 



from infinite plane to a bounded set, to a compact 
symmetric space, etc. But the elimination of para- 
dox destroys the essence of Carroll's problem. If we 
follow the paradox and try to give a meaning to 
"points are taken at random on an infinite plane" 
then we replace cj-additivity of the probability mea- 
sure by finite-additivity and come to the applied 
probability theory for finite-additive probabili- 
ties. Of course, this theory for abstract probability 
spaces would be too poor, and some additional ge- 
ometric and algebraic structures are necessary to 
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build rich enough theory. 

This is not just a beautiful geometrical problem, 
but rather an applied question about the proper def- 
inition of multiscale ensembles. We need such a def- 
inition to make any general statement about mul- 
tiscale systems, and briefly analyze lessons of Car- 
roll's problem in Sec. [31 

In this section we use some mathematics to define 
the multiscale ensembles with well separated con- 
stants. This is necessary background for the analy- 
sis of systems with limitation, and technical conse- 
quences are rather simple. We need only two proper- 
ties of a typical system from the multiscale ensemble 
with well separated constants: 

(i) Every two reaction rate constants fc, fc', are 
connected by the relation fc ^ fc' or fc ^ k' 
(with probability close to 1); 

(ii) The first property persists (with probability 
close to 1), if we delete two constants k, k' from 
the list of constants, and add a number fcfc' or 
a number k/k' to that list. 

If the reader can use these properties (when it is nec- 
essary) without additional clarification, it is possible 
to skip reading Sec. [3] and go directly to more ap- 
plied sections. In. Sec.[4]we study static and dynamic 
properties of linear multiscale reaction networks. An 
important instrument for that study is a hierarchy 
of auxiliary discrete dynamical system. Let Ai be 
nodes of the network ("components"), Ai Aj be 
edges (reactions) , and kji be the constants of these 
reactions (please pay attention to the inverse order 
of subscripts). A discrete dynamical system (f> is a 
map that maps any node Ai in a node ^^(i) • To con- 
struct a first auxiliary dynamical system for a given 
network we find for each Ai the maximal constant of 
reactions Ai — > Aj: k^(^i-^i > kji for all j, and (f)(i) = 
i if there are no reactions Ai Aj. Attractors in 
this discrete dynamical system are cycles and fixed 
points. 

The fast stage of relaxation of a complex reaction 
network could be described as mass transfer from 
nodes to correspondent attractors of auxiliary dy- 
namical system and mass distribution in the attrac- 
tors. After that, a slower process of mass redistribu- 
tion between attractors should play a more impor- 
tant role. To study the next stage of relaxation, we 
should glue cycles of the first auxiliary system (each 
cycle transforms into a point), define constants of 
the first derivative network on this new set of nodes, 
construct for this new network a (first) auxiliary dis- 
crete dynamical system, etc. The process terminates 
when we get a discrete dynamical system with one 



attractor. Then the inverse process of cycle restora- 
tion and cutting starts. As a result, we create an 
explicit description of the relaxation process in the 
reaction network, find estimates of eigenvalues and 
eigenvectors for the kinetic equation, and provide 
full analysis of steady states for systems with well 
separated constants. 

The problem of multiscale asymptotics of 
eigenval ues of non-selfadjoin t matrices wa s stud- 
i ed b y IVishik fc Ljusternikl (|l960l ) and iLidskiil 
Recently, some generalizations were ob- 
tained by idempotent (min-plus) algebra methods 
(jAkian. Bapat. fc Gauberd (|2004l) ). These methods 
provide a natural lan guage for discussion of som e 
multiscale problems ( Litvinov fc Maslov ( 2005h ). 
In the Vishik-Ljusternik-Lidskii theorem and its 
generalizations the asymptotics of eigenvalues and 
eigenvectors for the family of matrices Aij (e) = 
a^je^'^ + o(e'^'j ) is studied for e > 0, e 0. 

In the chemical reaction networks that we study, 
there is no small parameter e with a given dis- 
tribution of the orders e^^^ of the matrix nodes. 
Instead of these powers of e we have orderings of 
rate constants. Furthermore, the matrices of kinetic 
equations have some specific properties. The possi- 
bility to operate with the graph of reactions (cycles 
surgery) significantly helps in our constructions. 
Nevertheless, there exists some similarity between 
these problems and, even for general matrices, 
graphical representati on is useful. The language of 
idempotent algebra ( Litvinov fc Maslov ( 20051 )). 
as well as nonstandard analysis w ith infinitisemals 
( Albeverio. Fenstad, Hoegh-Krohn . fc L indstromI 
(Il986l )). can be used for description of the multi- 
scale reaction networks, but now we postpone this 
for later use. 

We summarize results of relaxation analysis and 
describe the algorithm of approximation of steady 
state and relaxation in Subsec. After that, sev- 
eral examples of networks are analyzed. In Sec. [5] 
we illustrate the analysis of dominant systems on 
a simple example, the reversible triangle of reac- 
tions: Ai ^ A2 A3 Ai. This simplest example 
became very popular for the lumping analysis case 
study after the well known work of iWei fc PrateJ 



(Il962l ). The most important mathematical proofs 
are presented in the appendices. 

In multiscale asymptotic analysis of reaction net- 
work we found several very attractive zero- one laws. 
First of all, components eigenvectors are close to 
or ±1. This law together with two other zero-one 
laws are discussed in Sec. (6] "Three zero-one laws 
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and nonequilibrium phase transitions in multiscale 
systems." 

A multiscale system where every two constants 
have very different orders of magnitude is, of course, 
an idealization. In parametric families of multiscale 
systems there could appear systems with several 
constants of the same order. Hence, it is necessary 
to study effects that appear due to a group of con- 
stants of the same order in a multiscale network. The 
system can have modular structure, with different 
time scales in different modules, but without sepa- 
ration of times inside modules. We discuss systems 
with modular structure in Sec.[7l The full theory of 
such systems is a challenge for future work, and here 
we study structure of one module. The elementary 
modules have to be solvable. That means that the 
kinetic equations could be solved in explicit analyti- 
cal form. We give the necessary and sufficient condi- 
tions for solvability of reaction networks. These con- 
ditions are presented constructively, by algorithm of 
analysis of the reaction graph. 

It is necessary to repeat our study for nonlinear 
networks. We discuss this problem and perspective 
of its solution in the concluding Sec.[8l Here we again 
use the ex perience summarized in the l UPAC Com- 
pendium ( Rate-controlling stepi (2007)) where the 
notion of controlling step is generalized onto non- 
linear elementary reaction by inclusion of some con- 
centration into "pseudo-first order rate constant" . 

2. Static and dynamic limitation in a linear 
chain and a simple catalytic cycle 

2.1. Linear chain 

A linear chain of reactions, Ai ^ A2 ^ ---An, 
with reaction rate constants ki (for Ai — > A^+i), 
gives the first example of limitation. Let the reac- 
tion rate constant kq be the smallest one. Then we 
expect the following behaviour of the reaction chain 
in time scale ~ l/fc^: all the components Ai, ...Aq^i 
transform fast into Ag, and all the components 
Aq+i, ...An-i transform fast into An, only two com- 
ponents, Aq and An are present (concentrations of 
other components are small) , and the whole dy- 
namics in this time scale can be represented by a 
single reaction Aq An with reaction rate con- 
stant kq . This picture becomes more exact when kq 
becomes smaller with respect to other constants. 

The kinetic equation for the linear chain is 



where Ci is concentration of Ai and fci_i ~ for i = 
1. The coefficient matrix K of this equations is very 
simple. It has nonzero elements only on the main 
diagonal, and one position below. The eigenvalues 
of K are —ki (i = 1, ...n ~ 1) and 0. The left and 
right eigenvectors for eigenvalue, f' and r°, are: 



/° = (1, !,...!), r" = (0,0,...0,l). 



(2) 



all coordinates of are equal to 1, the only nonzero 
coordinate of r° is and we represent vector- 
column in row. 

Below we use explicit form of K left and right 
eigenvectors. Let vector-column r* and vector-row 
r be right and left eigenvectors of K for eigenvalue 
— ki. For coordinates of these eigenvectors we use 
notation and Z*. Let us choose a normalization 
condition r\ — l\ = \. It is straightforward to check 
that rj = (j < i) and 1] = {j > i), r^^^^ = 
- ^i) (j > a-nd rj_^ kj^ilj/{kj^i - 
kj) (j < i), and 



n 



Vi+j-l 



j = l 



n 



ki^j 



(3) 



It is convenient to introduce formally fco = 0. Under 
selected normalization condition, the inner product 
of eigenvectors is: ZV-' = Sij, where Sij is the Kro- 
necker delta. 

If the rate constants are well separated (i.e., any 
two constants, ki, kj are connected by relation ki 3> 
kj or kj J 



ki—j ki 



1, if fci < ki-j] 



(4) 



0, iffe, >fe,_j . 

Hence, |Zi_„J « 1 or |Z,-_„| « 0. To demonstrate 
that also « 1 or |r-+,„| « 0, we shift nomina- 

tors in the product ^ on such a way: 



i+ra 



ki 



ki 



m— 1 

n 



+j 



— ki~-\Ci—\ k^Ci 



(1) 



Exactly as in (|4]) , each multiplier 
ki+j/{ki^j — ki) 
here is either almost 1 or almost 0, and 

is either almost or almost —1. In this zero-one 
asymptotics 
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1 



if > ki for all 7 = 1, ... m, else 

_1 ^ _1 

if fci+ j > ki for all = 1 , . . . m — 1 
and ki+ra < h, else r-_,_„ w 0. 



0; 



(5) 



In this asymptotic, only two coordinates of right 
eigenvector r* can have nonzero values, r\ = 1 and 
T\j^m « — 1 where m is the first such positive integer 
that i + m < n and fc,;+m < fc^. Such m always 
exists because fc„ = 0. For left eigenvector f , Z- « 
. . .U „ « 1 and li „ , « where i > and q is 
the first such positive integer that i — q — 1 > 
and ki-q-i < ki. It is possible that such q does 
not exist. In that case, all w 1 for j > 0. It 
is straightforward to check that in this asymptotic 

The simplest example gives the order fci 3> fc2 3> 
... > fc„_i: « 1 for J > 0, = 1, r]^, « -1 
and all other coordinates of eigenvectors are close to 
zero. For the inverse order, fci <C ^2 ^ •■• <Si fcn-i, 

= 1, = 1, « — 1 and all other coordinates of 
eigenvectors are close to zero. 

For less trivial example, let us find the asymptotic 
of left and right eigenvectors for a chain of reactions: 



where the upper index marks the order of rate con- 
stants: fc4 3> fcs 3> fc2 S> fca S> fci (fci is the rate 
constant of reaction Ai — s- ...). 

For left eigenvectors, rows Z', we have the following 
asymptotics: 

(1,0,0,0,0,0), (0,1,0,0,0,0), 

Z3«(0,1,1,0,0,0),Z4«(0,0,0,1,0,0), (6) 
(0,0,0,1,1,0). 

For right eigenvectors, columns r% we have the 
following asymptotics (we write vector-columns in 
rows): 

(1,0,0,0,0,-1), (0,1,-1,0,0,0), 
(0,0,1,0,0,-1), (0,0,0,1,-1,0), (7) 
(0,0,0,0,1,-1). 

The correspondent approximation to the general so- 
lution of the kinetic equations is: 

n-l 

c{t) = (/°c(0))r° + ^(rc(0))r^ exp(-fc,t), (8) 
1=1 



where c(0) is the initial concentration vector, and 
for left and right eigenvectors P and we use their 
zero-one asymptotic. 

Asymptotic formulas allow us to transform kinetic 
matrix X to a matrix with value of diagonal element 
could not be smaller than the value of any element 
from the correspondent column and row. 

Let us represent the kinetic matrix K in the ba- 
sis of approximations to eigenvectors ((T]) . The trans- 



formed matrix is Ki 



PKr^ = 0,1, ...5): 



K = 



K 



-fci 








fci ^ 


-fc2 








fc2 -fcs 








fca - 


-ki 








fc4 — fcs 








fcs 











-fc 


L 





fci 


-fc2 





fci 


fcs ~h 





-fc'3 fcs 


-fc4 





-fca fcs 


-fcs -^5 



(9) 



The transformed matrix has an important property 

\k,,\<m\n{\kul\K,,\}. 

The initial matrix K is diagonally dominant in 
columns, but its rows can include elements that 
are much bigger than the correspondent diagonal 
elements. 

We mention that a naive expectation Kij w 5ij 
is not realistic: some of the nondiagonal matrix ele- 
ments Kij are of the same order than To\n{Kii, Kjj}. 
This example demonstrates that a good approxima- 
tion to an eigenvector could be not an approximate 
eigenvector. If Ke = Ae and ||e — /|| is small then 
/ is an approximation to eigenvector e. If K f « A/ 
(i.e. \\K f — A/ 1 1 is small), then / is an approximate 
eigenvector for eigenvalue A. Our kinetic matrix K 
is very ill-conditioned. Hence, nobody can guaran- 
tee that an approximation to eigenvector is an ap- 
proximate eigenvector, or, inverse, an approximate 
eigenvector (a "quasimode" ) is an approximation to 
an eigenvector. 

The question is, what do we need for approxima- 
tion of the relaxation process ^ . The answer is ob- 
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vious: for approximation of general solution ([5]) with 
guaranteed accuracy we need approximation to the 
genuine eigenvectors ("modes") with the same ac- 
curacy. The zero-one asymptotic ([5|) gives this ap- 
proximation. Below we always find the modes ap- 
proximations and not quasimodes. 



2.2. General properties of a cycle 

The catalytic cycle is one of the most important 
substructures that we study in reaction networks. In 
the reduced form the catalytic cycle is a set of linear 
reactions: 



Ai A2 



.Ar. 



Ai. 



Reduced form means that in reality some of these 
reaction are not monomolecular and include some 
other components (not from the list Ai, . . . A„). But 
in the study of the isolated cycle dynamics, concen- 
trations of these components are taken as constant 
and are included into kinetic constants of the cycle 
linear reactions. 

For the constant of elementary reaction Ai we 
use the simplified notation ki because the product 
of this elementary reaction is known, it is ^i+i for 
i < n and Ai for i — n. The elementary reaction 
rate is Wi = kiCi, where Ci is the concentration of 
Ai. The kinetic equation is: 

Ci = Wi-i - Wi, (10) 

where by definition wq = w„. In the stationary state 
(ci = 0), all the Wi are equal: Wi = w. This common 
rate w we call the cycle stationary rate, and 



w 



1 



W 
ki 



(11) 



where b = Ci is the conserved quantity for reac- 
tions in constant volume (for general case of chem- 
ical kinetic equations see elsewhere, for example, 
the bo ok bv lYablonskii. Bvkov. Gorban. &: Elokhin 
(|l99l[ )). The stationary rate w (fTTj) is a product of 
the arithmetic mean of concentrations, b/n, and the 
harmonic mean of constants (inverse mean of inverse 

h). 



2.3. Static limitation in a cycle 



If one of the constants, /cmin, is much smaller than 




(12) 



w = knb, (13) 



or simply in linear approximation 

\ i<n / 

where we should keep the first-order terms in c„ in 
order not to violate the conservation law. 

The simplest zero order approximation for the 
steady state gives 

c„ = 6, c^^O{iy^ n). (14) 

This is trivial: all the concentration is collected at 
the starting point of the "narrow place" , but may be 
useful as an origin point for various approximation 
procedures. 

So, the stationary rate of a cycle is determined 
by the smallest constant, fc„iin, if fcmin is sufficiently 
small: 

w^k^,^b if (15) 

In that case we say that the cycle has a limiting step 
with constant fcmin- 



2.4. Dynamical limitation in a cycle 

If kn/ki is small for all i < n, then the kinetic be- 
haviour of the cycle is extremely simple: the coeffi- 
cients matrix on the right hand side of kinetic equa- 
tion pn]) has one simple zero eigenvalue that corre- 
sponds to the conservation law ^ c,; — b and n — 1 
nonzero eigenvalues 



Aj = -fcj + Si {i < n). 



(16) 



others (let it be k„ 



kn), then 



where 6i ^ when J2i<n 17 ^ 0. 

It is easy to demonstrate p^ : let us exclude the 
conservation law (the zero eigenvalue) ^Ci = b and 
use independent coordinates q {i — 1, . . .n — 1); 
Cn = b — J2i<n ^i- these coordinates the kinetic 
equation (fTO|) has the form 

c = Kc- knAc + knbe^ (17) 

where c is the vector-column with components Ci 
(i < n), K is the lower triangle matrix with nonzero 
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elements only in two diagonals: {K)ii — —ki (i = 
1, . . . n - 1), {K)i+i^ i = ki {i = 1,. . .n - 2) (this 
is the kinetic matrix for the linear chain of n — 1 
reactions Ai ^ A2 ^ •••^n); A is the matrix with 
nonzero elements only in the first row: = 1, 

is the first basis vector (e} = l,ej = Oforl <i<n). 
After that, eq. (jl6p follows simply from continuous 
dependence of spectra on matrix. 

The relaxation time of a stable linear system (I17p 
is, by definition, 



T = [min{i?e(— Ai) | i = 1, . . . 71 — 1}] 



-1 



For small fc„, 

TKil/kr, kr = Taiii{ki\i = I, . . .n — 1}. (18) 
In other words, k^ is the second slowest rate con- 

2.5. Relaxation equation for a cycle rate 

A definition of the cycle rate is clear for steady 
states because stationary rates of all elementary re- 
actions in cycle coincide. There is no common def- 
inition of the cycle rate for nonstationary regimes. 
In practice, one of steps is the step of product re- 
lease (the "final" step of the catalytic transforma- 
tion), and we can consider its rate as the rate of the 
cycle. Formally, we can take any step and study re- 
laxation of its rate to the common stationary rate. 
The single relaxation time approximation gives for 
rate Wi of any step: 



= kr{k^inb - Wi); 

Wi{t) = fcminfc + e-^^\w^{Q) - k^ij)), 



(19) 



where kmi-a is the limiting (the minimal) rate con- 
stant of the cycle, k-r is the second in order rate con- 
stant of the cycle. 

So, for catalytic cycles with the limiting constant 
/cniinj the relaxation time is also determined by one 
constant, but another one. This is kr, the second 
in order rate constant. It should be stressed that 
the only smallness condition is required, /cmin should 
be much smaller than other constants. The second 
constant, k^- should be just smaller than others (and 
bigger than fcmin), but there is no ^ condition for 
k-r required. 

One of the methods for measurement of chemi- 
cal reaction c onstants is the relaxation spectroscopy 
(Eigen ( 19721 )). Relaxation of a system after an im- 
pact gives us a relaxation time or even a spectrum of 
relaxation times. For catalytic cycle with limitation. 



the relaxation experiment gives us the second con- 
stant kr, while the measurement of stationary rate 
gives the smallest constant, k-cai-a- This simple re- 
mark may be important for relaxation spectroscopy 
of open system. 

2.6. Ensembles of cycles and robustness of 
stationary rate and relaxation time 

Let us consider a catalytic cycle with random rate 
constants. For a given sample constants /ci, . . . fc„ 
the zth order statistics is equal its zth-smallest value. 
We are interested in the first order (the minimal) 
and the second order statistics. 

For independent identically distributed constants 
the variance of fcmin = minjfci, . . . is signifi- 
cantly smaller then the variance of each fc^, Var(A;). 
The same is true for statistic of every order. For 
many important distributions (for example, for 
uniform distribution), the variance of ith order 
statistic is of order ~ \ai{k)/n^ . For big n it goes 
to zero faster than variance of the mean that is of 
order ~ Var(fc)/n. To illustrate this, let us consider 
n constants distributed in interval [a, 6]. For each 
set of constants, /ci, . . . fc„ we introduce "symmetric 
coordinates" s^: first, we order the constants, a < 



< < 



a, 



"0 — ■"ij-i-i 
Transformation 



^ < b, then calculate sq — ki-^ 
(j = 1, ... n - 1), Sn ^ b - h^. 
(fci,...fc„) 1-^ (so,...s„) maps 
a cube [a, 6]" onto n-dimensional simplex A„ = 
{(so, ■ • ■ Sn) I X^i Si = b — a} and uniform distribu- 
tion on a cube transforms into uniform distribution 
on a simplex. 

For large n, almost all volume of the simplex is 
concentrated in a small neighborhood of its center 
and this effect is an example of measure concen- 
tration effects that pla y important rol e in modern 
geometry and analysis ( Gromov ( 1999t )). All Si are 
identically distributed, and for normalized variable 
s = Si/{b — a) the first moments are: E(s) = l/{n + 



1) = l/n + o{l/n), 
2/7i2 + o(l/7i2), 



E(s^ 



2/[(n+ l)(n + 2)] 



Var(s) =E(s2) - (E(s))2 



(n + l)2(n-K2) 



7^2 



Hence, for example, Var(fcinin) — [b - 
o{l/n'^). The standard deviation of fcmin goes to 
zero as 1/n when n increases. This is much faster 
than l/\/n prescribed to the deviation of the mean 
value of independent observation (the "law of er- 



9 



rors"). The same asymptotic 1/n is true for the 
standard deviation of the second constant also. 
These parameters fluctuate much less than individ- 
ual constants, and even less than mean constant 
(for more examples with applicati ons to statistic al 
physics we address to the paper bv lGorban ( 2006f l). 

It is impossible to use this observation for cy- 
cles with limitation directly, because the inequal- 
ity of limitation (llSp is not true for uniform distri- 
bution. According to this inequality, ratios ki/k-arm 
should be sufficiently small (if ki ^ fcmin). To pro- 
vide this inequality we need to use at least the log- 
uniform distribution: fc; = exp and are inde- 
pendent variables uniformly distributed in interval 
[a, (3\ with sufficiently big (/3 — a)/n. 

One can interpret the log-uniform distribution 
through the Arrhenius law: k — 74exp(— AG/fcT), 
where AG is the change of the Gibbs free energy 
inreaction (it includes both energetic and entropic 
terms: AG = Ai/ - TAS", where Ai/ is enthalpy 
change and AS" is entropy change in reaction, T is 
temperature) . The log- uniform distribution of k cor- 
responds to the uniform distribution of AG. 

For log- uniform distribution of constants fci , . . . 
kn, if the interval of distribution is sufficiently big 
(i.e. {(3 — a)/n ^ 1), then the cycle with these 
constants has the limiting step with probability 
close to one. More precisely we can show that for 
any two constants ki,kj the probability l?[ki/kj > 
r or kj/ki > r] = (1 — \og{r)/{(3 — a))^ approaches 
one for any fixed r > 1 when (3 — a — > oo. Re- 
laxation time of this cycle is determined by the 
second constant kr (also with probability close to 
one). Standard deviations of /cmin and kr are much 
smaller than standard deviation of single constant 
ki and even smaller than standard deviation of mean 
constant ki/n. This effect of stationary rate and 
relaxation time robustness seems to be important 
for understanding robustness in biochemical net- 
works: behaviour of the entire system is much more 
stable than the parameters of its parts; even for 
large fluctuations of parameters, the system does 
not change significantly the stationary rate (statics) 
and the relaxation time (dynamics) . 



2.7. Systems with well separated constants and 
monotone relaxation 

The log-uniform identical distribution of indepen- 
dent constants ki, . . . kn with sufficiently big inter- 
val of distribution {{(3 — a)/n^ 1) gives us the first 



example of ensembles with well separated constants: 
any two constants are connected by relation 3> or 
<^ with probability close to one. Such systems (not 
only cycles, but much more complex networks too) 
could be studied analytically "up to the end" . 

Some of their properties are simpler than for gen- 
eral networks. For example, the damping oscillations 
are impossible, i.e. the eigenvalues of kinetic matrix 
are real (with probability close to one) . If constants 
are not separated, damped oscillations could exist, 
for example, if all constants of the cycle are equal, 
ki = k2 = ■ ■ ■ = kn = k, then (1 + A/fc)" = 1 and 
Am = /c(exp(27rim/n) — 1) (m = 1, ... n — 1), the 
case m = corresponds to the linear conservation 
law. Relaxation time of this cycle may be relatively 
big: T = i(l-cos(27r/n))-i - n^/{27Tk) (for big n). 

The catalytic cycle without limitation can have 
relaxation time much bigger then l/Zcmin, where 
/cniin is the minimal reaction rate constant. For ex- 
ample, if all k are equal, then for ti = 11 we get 
T « 20/fc. In more detail the possible relations be- 
tw een r and the slowest consta nt were discussed 
bv lYablonskii &: Cheresi j (1984). In that paper, a 
variety of cases with different relationships between 
the steady-state reaction rate and relaxation was 
presented. 

For catalytic cycle, if a matrix K~knA (fT7|) has a 
pair of complex eigenvalues with nonzero imaginary 
part, then for some g € [0, 1] the matrix K — gknA 
has a degenerate eigenvalue (we use a simple con- 
tinuity argument). With probability close to one, 
^min ^ ~ for ^''^Y two ki , kj that are not 
minimal. Hence, the fcmin-small perturbation can- 
not transform matrix K with eigenvalues ki p6p 
and given structure into a matrix with a degenerate 
eigenvalue. For proof of this statement it is sufficient 
to refer to diagonal dominance of K (the absolute 
value of each diagonal element is greater than the 
sum of the absolute values of the other elements in 
its column) and classical inequalities. 

The matrix elements of A in the eigenbasis of 
K are {A)ij — VAr^ . From obtained estimates for 
eigenvectors we get |(^)jj| ^ 1 (with probability 
close to one). This estimate does not depend on 
values of kinetic constants. Now, we can apply the 
G ershgorin theorem (see , for example, the review 
of 'Marcus & Mine (1992^ and for more details the 
book of Varga (2004)) to the matrix K — knA in the 
eigenbasis of K: the characteristic roots of if — knA 
belong to discs |z -I- < knRi{A), where Ri{A) = 
J2j I {A)ij I . If the discs do not intersect, then each of 
them contains one and only one characteristic num- 
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ber. For ensembles with well separated constants 
these discs do not intersect (with probability close 
to one). Complex conjugate eigenvalues could not 
belong to different discs. In this case, the eigenval- 
ues are real - there exist no damped oscillations. 

2.8. Limitation by two steps with comparable 
constants 

If we consider one-parametric families of systems, 
then appearance of systems with two comparable 
constants may be unavoidable. Let us imagine a con- 
tinuous path ki{s) (s e [0, 1], s is a parameter along 
the path) in the space of systems, which goes from 
one system with well separated constants (s = 0) to 
another such system (s = 1). On this path ki{s) such 
a point s that fci(s) = kj{s) may exist, and this ex- 
istence may be stable, that is, such a point persists 
under continuous perturbations. This means that on 
a path there may be points where not all the con- 
stants are well separated, and trajectories of some 
constants may intersect. 

For catalytic cycle, we are interested in the follow- 
ing intersection only: fcmin and the second constant 
are of the same order, and are much smaller than 
other constants. Let these constants be kj and ki, 
j ^ I. The limitation condition is 
1 



(20) 



The steady state reaction rate and relaxation time 
are determined by these two constants. In that case 
their effects are coupled. For the steady state we get 
in first order approximation instead of (|13p : 



Cl 




kj + ki 



(21) 

Elementary analysis shows that under the limitation 
condition (PO)) the relaxation time is 

kj + ki ^ ^ 

The single relaxation time approximation for all 
elementary reaction rates in a cycle with two limit- 
ing reactions is 



Wi = kjkib - {kj + kijWi] 

(23) 

The catalytic cycle with two limiting reactions has 
the same stationary rate w (j21[) and relaxation time 
([22]) as a reversible reaction A ^ B with fc"*" = kj, 
k^ ~ ki. 

In two-parametric families three constants can 
meet. If three smallest constants kj , ki , km have com- 
parable values and are much smaller than others, 
then static and dynamic properties would be deter- 
mined by these three constants. Stationary rate w 
and dynamic of relaxation for the whole cycle would 
be the same as for 3-reaction cycle A ^ B ^ C ^ 
A with constants kj , ki , km ■ The damped oscillation 
here are possible, for example, if kj = ki = km = k, 
then there are complex eigenvalues A = A;(— |±i-^). 
Therefore, if a cycle manifests damped oscillation, 
then at least three slowest constants are of the same 
order. The same is true, of course, for more general 
reaction networks. 

In iV-parametric families of systems N+1 smallest 
constants can meet, and near such a "meeting point" 
a slow auxiliary cycle of iV -I- 1 reactions determines 
behaviour of the entire cycle. 



2.9. Irreversible cycle with one inverse reaction 

In this subsection, we represent a simple example 
that gives the key to most of subsequent construc- 
tions of "cycles surgery" . Let us add an inverse re- 
action to the irreversible cycle: Ai ... Ai ^ 
Ai+i ... An A\. We use the previous no- 
tation k\,...kn for the cycle reactions, and k~ for 
the inverse reaction Ai ^ A^+i. For well-separated 
constants, influence of k'^ on the whole reaction is 
determined by relations of three constants: fc^, k~ 
and /ci+i. First of all, if k'^ <^ ki^i then in the main 
order there is no such influence, and dynamic of the 
cycle is the same as for completely irreversible cycle. 

If the opposite inequality is true, k^ ^ fci+i, 
then equilibration between Ai and Ai^i gives 
kiCiX w k^Ci-^-i. If we introduce a lumped compo- 
nent Aj with concentration cj = c; -I- Q+i, then 
Cl w k^cj/ {ki + kl) and a+i « hcj/ih + k~). Us- 
ing this component instead of the pair Ai,Ai+i we 
can consider an irreversible cycle with n — 1 com- 
ponents and n reactions Ai — > ... — > Ai-i A\ ^ 
Ai^2 An — > ^1. To estimate the reaction 

rate constant k} for a new reaction, A\ Ai^2-, 
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let us mention that the correspondent reaction rate 
should be fci+iCi+i ~ ki^ikic] / (ki + k^). Hence, 



Ai 9 — 



- (ki + fcj + ki+i) 



ki/ {ki 



For systems with well separated constants this ex- 
pression can be simplified: if ki 3> k^ then kj w fc^+i 
and if h <C k~ then kj « ki+ih/k^ . The first case, 
ki ^ k~ is limitation in the small cycle (of length 
two) Ai <-> Ai+i by the inverse reaction Ai ^ Ai+i. 
The second case, fc^ <C fe~, means the the direct re- 
action is the limiting step in this small cycle. 

In order to estimate eigenvectors, we can, after 
identification of the limiting step in the small cycle, 
delete this step and reattach the outgoing reaction 
to the beginning of this step. For the first case, ki ^ 
k~ , we get the irreversible cycle, Ai ... Ai ^ 
Ai+i ... An Ai, with the same reaction rate 
constants. For the second case, ki <C k~ we get a new 
system of reactions: a shortened cycle Ai ... 
Ai Ai+2 ^ ... An ^ Ai and an "appendix" 
Ai^i Ai. For the new elementary reaction Ai 
Ai+2 the reaction rate constant is k} « ki^iki/k^ . 
All other elementary reactions have the same rate 
constants, as they have in the initial system. After 
deletion of the limiting step from the "big cycle" 
Ai ... Ai ^ Ai+2 ■■■ An Ai, we 
get an acyclic system that approximate relaxation 
of the initial system. 

So, infiuence of a single inverse reaction on the 
irreversible catalytic cycle with well-separated con- 
stants is determined by relations of three constants: 
ki, k~ and fcj+i. If k~ is much smaller than at least 
one of ki, ki^i, then there is no influence in the main 
order. If k~ ^ ki and k~ ^ fc^+i then the relax- 
ation of the initial cycle can be approximated by re- 
laxation of the auxiliary acyclic system. 

Asymptotic equivalence (for k~ ^ ki, fci+i) of 
the reaction network Ai ^ Ai^i Ai^2 with rate 
constants ki, k~ and ki+i to the reaction network 
Ai+i Ai ^ Ai+2 with rate constants k~ (for the 
reaction A^+i Ai) and ki+iki/k~ (for the reac- 
tion Ai — !■ Ai_|_2) is simple, but slightly surprising 
fact. The kinetic matrix for the first network in co- 
ordinates Ci, Ci+i, is 



K 



-h k~ 
h -{kl 
fc.+i 



Ai = -fc.+ifc.(l + o(l))/fc,-, A2 = -kr{l + o(l)), 
where o(l) <C 1. Right eigenvector for zero eigen- 
value is (0, 0, 1) (we write vector columns in rows). 
For Ai the eigenvector is = (1, 0, —1) + o(l), and 
for A2 it is — (1,-1, 0) +o(l). For the linear chain 
of reactions, Ai^i Ai —^ Ai^2, with rate con- 
stants k~ and ki^iki/k~ eigenvalues are —k^ and 
-~ki^iki/k~ . These values approximate eigenvalues 
of the initial system with small relative error. The 
linear chain has the same zero-one asymptotic of the 
correspondent eigenvectors. 

This construction, a small cycle inside a big sys- 
tem, a quasi steady state in the small cycle, and 
deletion of the limiting step with reattaching of reac- 
tions (see Fig. [l]below) appears in this paper many 
times in much general settings. The uniform esti- 
mates that we need for approximation of eigenvalues 
and eigenvectors by these procedures are proven in 
Appendices. 

3. Multiscale ensembles and finite additive 
distributions 

3.1. Ensembles with well separated constants, 
formal approach 

In previous section, ensembles with well separated 
constants appear. We represented them by a log- 
uniform distribution in a sufficiently big interval 
logfc e but we were not interested in most 

of probability distribution properties, and did not 
use them. The only property we really used is: if 
ki > kj, then ki/kj ^ 1 (with probability close to 
one). It means that we can assume that ki/kj ^ a 
for any preassigned value of a that does not depend 
on k values. One can interpret this property as an 
asymptotic one for a — > — 00, /? — > 00. 

That property allows us to simplify algebraic for- 
mulas. For example, ki + kj could be substituted by 
maxjfci, kj} (with small relative error), or 



aki + bkj 



( a/c, 
\b/d. 



if ki > k 
if ki <^ ki 



J' 



The eigenvalues are and 



for nonzero a, b, c, d (see, for example, 

Of course, some ambiguity can be introduced, for 
example, what is it, (fci -I- ^2) — fci, if fci ^ fc2? If we 
first simplify the expression in brackets, it is zero, 
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but if we open brackets without simplification, it is 
k2 ■ This is a standard difficulty in use of relative er- 
rors for round-off. If we estimate the error in the final 
answer, and then simplify, we shall avoid this diffi- 
culty. Use of o and O symbols also helps to control 
the error qualitatively: if fci 3> fc2, then we can write 
(fci + fca) = fci(l + o(l)), and fci(l + o(l)) - fci = 
fcio(l). The last expression is neither zero, nor abso- 
lutely small - it is just relatively small with respect 
to ki. 

The formal approach is: for any ordering of rate 
constants, we use relations ^ and <C and assume 
that ki/kj ^ a for any preassigned value of a that 
does not depend on k values. This approach allows 
us to perform asymptotic analysis of reaction net- 
works. A special version of this approach consists of 
group ordering: constants are separated on several 
groups, inside groups they are comparable, and be- 
tween groups there are relations ^ or <C. An exam- 
ple of such group ordering was discussed at the end 
of previous section (several limiting constants in a 
cycle). 



3.2. Probability approach: finite additive measures 

The asymptotic analysis of multiscale systems for 
log-uniform distribution of independent constants 
on an interval log k e [a, (3] (—a, (3 — *■ oo) is possi- 
ble, but parameters a, [3 do not present in any an- 
swer, they just should be sufficiently big. A natural 
question arises, what is the limit? It is a log-uniform 
distribution on a line, or, for n independent identi- 
cally distributed constants, a log-uniform distribu- 
tion on M"). 

It is well known that the uniform distribution on 
M" is impossible: if a cube has positive probabil- 
ity e > (i.e. the distribution has positive density) 
then the union of > 1/e such disjoint cubes has 
probability bigger than 1 (here we use the finite- 
additivity of probability). This is impossible. But if 
that cube has probability zero, then the whole space 
has also zero probability, because it can be covered 
by countable family of the cube translation. Hence, 
translation invariance and cr-additivity (countable 
additivity) are in contradiction (if we have no doubt 
about probability normalization). 

Nevertheless, there exists finite-additive proba- 
bility which is invariant with respect to Euclidean 
group E[n) (generated by rotations and transla- 
tions). Its values are densities of sets. 

Let A be the Lebesgue measure in M" and D C M" 



be a Lebesgue measurable subset. Density of D is 
the limit (if it exists) : 

\{Dr\W;:) 



p{D) = lim 



(24) 



a(b;0 ' 

where B" is a ball with radius r and centre at origin. 
Density of M" is 1, density of every half-space is 1/2, 
density of bounded set is zero, density of a cone is its 
solid angle (measured as a sphere surface fractional 
area). Density ([24|) and translation and rotational 
invariant. It is finite-additive: if densities p{D) and 
p{H) ^ exist and DnH = then p{DUH) exists 
and p{D UH)= p{D) + p{H). 

Every polyhedron has a density. A polyhedron 
could be defined as the union of a finite number of 
convex polyhedra. A convex polyhedron is the inter- 
section of a finite number of half-spaces. It may be 
bounded or unbounded. The family of polyhedra is 
closed with respect to union, intersection and sub- 
traction of sets. For our goals, polyhedra form suffi- 
ciently rich class. It is important that in definition of 
polyhedron finite intersections and unions are used. 
If one uses countable unions, he gets too many sets 
including all open sets, because open convex poly- 
hedra (or just cubes with rational vertices) form a 
basis of standard topology. 

Of course, not every measurable set has density. If 
it is ne c essary , we can use the Hahn-Banach theorem 
RudinI (|l991[ ) and study extensions pEx of p with the 
following property: 



p{D) < PEx(-D) < piD), 



where 



p{D) — liminf 

— r — >oo 

p(Z3) = lim sup 



A(i:>nBp) 
A(£inB;f) 



A(B^) ■ 

Functionals p{D) and 75(13) are defined for all mea- 
surable D. We should stress that such extensions 
are not unique. Extension of density (j24p using the 
Hahn-Banach theorem for picking up a random in- 
teger was used in a very recent work bv lAdamaszek 
(l2006h . 



One of the most important concepts of any prob- 
ability theory is the conditional probability. In the 
density-based approach we can introduce the con- 
ditional density. If densities p{D) and p{H) ex- 
ist, p{H) and the following limit p{D\H) exists, 
then we call it conditional density: 

A(Dni? ni;?) 



p{D\H) = lim 



A(i7 nip) 



(25) 
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For polyhedra the situation is similar to usual prob- 
ability theory: densities p{D) and p{H) always exist 
and if p{H) ^ then conditional density exists too. 
For general measurable sets the situation is not so 
simple, and existence of piD) and p(-ff) 7^ does 
not guarantee existence of p{D\H). 

On a line, convex polyhedra are just intervals, fi- 
nite or infinite. The probability defined on polyhe- 
dra is: for finite intervals and their finite unions it is 
zero, for half-lines x>Q;ora:<Q!itisl/2, and for 
the whole line R the probability is 1. If one takes a set 
of positive probability and adds or subtracts a zero- 
probability set, the probability does not change. 

If independent random variables x and y are uni- 
formly distributed on a line, then their linear combi- 
nation z = ax -\- (3y is also uniformly distributed on 
a line. (Indeed, vector [x, y) is uniformly distributed 
on a plane (by definition), a set 2: > 7 is a half-plane, 
the correspondent probability is 1/2.) This is a sim- 
ple, but useful stability property. We shall use this 
result in the following form. If independent random 
variables /ci , . . . fc„ are log- uniformly distributed on 
a line, then the monomial J^"^]^ fc"' for real is 
also log-uniformly distributed on a line, if some of 
a., ^ 0. 



are absolutely no reasons for uniformity of C dis- 
tribution. And it is more important that the abso- 
lutely standard reasoning for independently chosen 
points gives another answer than could be found on 
the base of joint distribution. Why these approaches 
are in disagreement now? Because there is no classi- 
cal Fubini theorem for our finite-additive probabil- 
ities, and we cannot easily transfer from a multiple 
integral to a repeated one. 

There exists a much simpler example. Let x and 
y be independent positive real number. This means 
that vector (x, y) is uniformly distributed in the first 
quadrant. What is probability that x > yl Following 
the definition of probability based on the density 
of sets, we take the correspondent angle and find 
immediately that this probability is 1/2. This meets 
our intuition well. But let us take the first number x 
and look for possible values of y. The result: for given 
X the second number y is uniformly distributed on 
[0,00), and only a finite interval [0,x] corresponds 
to X > y. For the infinite rest we have x < y. Hence, 
X < y with probability 1. This is nonsense because 
of symmetry. So, for our finite-additive measure we 
cannot use repeated integrals (or, may be, should 
use them in a very peculiar manner). 



3.3. Carroll's obtuse problem and paradoxes of 
conditioning 



Lew is Carroll's Pillow Problem #58 (Carroll 
(jl958[ )): "Three points are taken at random on an 
infinite plane. Find the chance of their being the 
vertices of an obtuse-angled triangle." 

A random triangle on an infinite plane is presented 
by a point equidistributed in R^. Due to the den- 
sity - based definition, we should take and calculate 
the density of the set of obtuse-angled triangles in 
R®. This is equivalent to the problem: find a fraction 
of the sphere §^ C R® that corresponds to obtuse- 
angled triangles. Just integrate... . But there re- 
mains a problem. Vertices of triangle are indepen- 
dent. Let us use the standard logic for discussion of 
independent trials: we take the first point A at ran- 
dom, then the second point B, and then the third 
point C. Let us draw the first side AB. Immediately 
we find that for almost all positions of the the third 
point C the triangle is obtuse-angled (Guy (1993)). 
L. Carroll proposed to take another condition: let 
AB be the longest side and let C be uniformly dis- 
tributed in the allowed area. The answer then is easy 
— just a ratio of areas of two simple figures. But there 



3.4. Law of total probability and orderings 

For polyhedra, there appear no conditioning prob- 
lems. The law of total probabilities holds: if R" = 
U™ iJj are polyhedra, p{H,) > 0, p{H, nHj) = 
for i ^ j, and D C R" is a polyhedron, then 

rn rn 

p{D) = ^ p(I? n H,) = PiDmPm. (26) 

1=1 1=1 

Our basic example of multiscale ensemble is log- 
uniform distribution of reaction constants in R" 
(log ki are independent and uniformly distributed on 
the line). For every ordering kj-^ > kj^ > . . . > kj^ 
a polyhedral cone Hjij2...jn '^^ defined. These 

cones have equal probabilities p{Hj^j2...j^) = 1/nl 
and probability of intersection of cones for different 
orderings is zero. Hence, we can apply the law of to- 
tal probability (|26p . This means that we can study 
every event D conditionally, for different orderings, 
and then combine the results of these studies in the 
final answer ([26|) . 

For example, if we study a simple cycle then for- 
mula (jl3p for steady state is valid with any given 
accuracy with probability one for any ordering with 
the given minimal element fc„. 
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For cycle with given ordering of constants we can 
find zero-one approximation of left and right eigen- 
vectors ([5]). This approximation is valid with any 
given accuracy for this ordering with probability 
one. 

If we consider sufficiently wide log-uniform distri- 
bution of constants on a bounded interval instead of 
the infinite axis then these statements are true with 
probability close to 1. 

For general system that we study below the situ- 
ation is slightly more complicated: new terms, aux- 
iliary reactions with monomial rate constants k^^ = 
kl' could appear with integer (but not necessary 
positive) and we should include these in or- 
dering. It follows from stability property that these 
monomials are log-uniform distributed on infinite 
interval, if ki are. Therefore the situation seems to 
be similar to ordering of constants, but there is a sig- 
nificant difference: monomials are not independent, 
they depend on ki with q ^ 0. 

Happily, in the forthcoming analysis when we in- 
clude auxiliary reactions with constant k^, we al- 
ways exclude at least one of reactions with rate con- 
stant ki and q ^ 0. Hence, for we always can use the 
following statement (for the new list of constants, 
or for the old one): if kj^ > %2 > • • • > kj^ then 
kj-^ ^ fcj2 ^ . . . ^ kj^ , where a ^ b for positive 
a, b means: for any given e > the inequality ea > b 
holds with unite probability. 

If we use sufficiently wide but finite log-uniform 
distribution then e could not be arbitrarily small 
(this depends on the interval with), and probability 
is not unite but close to one. For given e > prob- 
ability tends to one when the interval width goes to 
infinity. It is important that we use only finite num- 
ber of auxiliary reactions with monomial constants, 
and this number is bounded from above for given 
number of elementary reactions. For completeness, 
we should mention here general algebraic theory of 
orderi ngs that i s nece s sary in more sophistica ted 
cases (|Robbianol (|l985h : ICreuel fc Pfistei] ('2002^ 'I. 

Finally, what is a most general form of multiscale 
ensemble? Our basic example of is the log-uniform 
distribution of reaction constants: logfc^ {i = l,n) 
are independent and uniformly distributed on the 
line. Of course, the main candidate to the general 
definition of multiscale ensemble is a distribution 
which is absolutely continuous with respect to that 
log-uniform distribution. If fi and v are measures 
on the same algebra of sets, then /i is said to be 
absolutely continuous with respect to v if fJ.{A) — 
for every set A for which i'(A) = 0. It is written as 



"/^ < i^" : 



HA) = =^ ^^{A) = 0) , 



For example, let us take an event (an infinite poly- 
hedron in logarithmic coordinates) U G i?" with 
nonzero probability in this distribution. The condi- 
tional probability distribution under condition U is 
an multiscale ensemble. We have such an ensemble 
for each polyhedral cone U with non-empty interior 
in the space of logarithms of constants. 

4. Relaxation of multiscale networks and 
hierarchy of auxiUary discrete dynamical 
systems 

4.1. Definitions, notations and auxiliary results 

4.1.1. Notations 

In this Sec, we consider a general network of lin- 
ear (monomolecular) reactions. This network is rep- 
resented as a directed graph (digraph): vertices cor- 
respond to components Ai, edges correspond to re- 
actions Ai Aj with kinetic constants kji > 0. For 
each vertex, Ai , a positive real variable Ci (concen- 
tration) is defined. A basis vector corresponds to 
Ai with components = Sij , where 6ij is the Kro- 
necker delta. The kinetic equation for the system is 

^^(fcyCj — fcjiCi), (27) 



dt 



or in vector form: c = Kc. 

To write another form of (f27|) we use stoichiomet- 
ric vectors: for a reaction Ai Aj the stoichiomet- 
ric vector is a vector in concentration space with 
ith coordinate —1, jth coordinate 1, and zero other 
coordinates. The reaction rate Wj 
netic equation has the form 



kjiCi- The ki- 



dc \ - 



(28) 



2J 



where c is the concentration vector. One more form 
of ((27| describes directly dynamics of reaction rates: 



dt 



I. ^ 
dt 



{wii - wii). (29) 



It is necessary to mention that, in general, system 
([29|) is not equivalent to (|28p. because there are ad- 



ditional connections between variables Wji. If there 



exists at least one Ai with two different outgoing 
reactions, Ai Aj and Ai Ai {j ^ I), then 
Wji/wii = kji/kii. If the reaction network generates 
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a discrete dynamical system Ai — > Aj on the set of 
Ai (see below), then the variables Wji are indepen- 
dent, and (j29p gives equivalent representation of ki- 
netics. 

For analysis of kinetic systems, linear conserva- 
tion laws and positively invariant polyhedra are im- 
portant. A linear conservation law is a linear func- 
tion defined on the concentrations 6(c) = '^i^iCi, 
whose value is preserved by the dynamics ([27|) . The 
conservation laws coefficient vectors hi are left eigen- 
vectors of the matrix K corresponding to the zero 
eigenvalue. The set of all the conservation laws forms 
the left kernel of the matrix K. Equation (p7)) al- 
ways has a linear conservation law: b^{c) = Ci = 
const. If there is no other independent linear conser- 
vation law, then the system is weakly ergodic. 

A set E is positively invariant with respect to ki- 
netic equations ^7\ . if any solution c{t) that starts 
in E at time to {c{to) £ E) belongs to E for t > to 
{c{t) £ E ii t > to). It is straightforward to check 
that the standard simplex E = {c | q > 0, q — 
1} is positively invariant set for kinetic equation 
(|27p : just to check that if Ci = for some i, and all 
Cj > then > 0. This simple fact immediately 
implies the following properties of K: 

- All eigenvalues A of X have non-positive real 
parts, ReX < 0, because solutions cannot leave E 
in positive time; 

- If ReX = then A = 0, because intersection of E 
with any plain is a polygon, and a polygon can- 
not be invariant with respect of rotations to suf- 
ficiently small angles; 

- The Jordan cell of K that corresponds to zero 
eigenvalue is diagonal - because all solutions 
should be bounded in E for positive time. 

- The shift in time operator exTp{Kt) is a contrac- 
tion in the li norm for t > 0: for positive t and 
any two solutions of ((27|) c{t),c'{t) € E 

J2hit)~c'M\<J2\^^iO)-c'M\■ 

i i 

Two vertices are called adjacent if they share a 
common edge. A path is a sequence of adjacent ver- 
tices. A graph is connected if any two of its vertices 
are linked by a path. A maximal connected sub- 
graph of graph G is called a connected component of 
G. Every graph can be decomposed into connected 
components. 

A directed path is a sequence of adjacent edges 
where each step goes in direction of an edge. A ver- 
tex A is reachable by a vertex B, if there exists an 
oriented path from B to A. 



A nonempty set V of graph vertexes forms a 
sink, if there are no oriented edges from Ai £ V to 
any Aj (jz. V. For example, in the reaction graph 
Ai ^ A2 ^ A3 the one-vertex sets {Ai} and {A3} 
are sinks. A sink is minimal if it does not contain 
a strictly smaller sink. In the previous example, 
{Ai}, {A3} are minimal sinks. Minimal sinks are 
also called ergodic components. 

A digraph is strongly connected, if every vertex 
A is reachable by any other vertex B. Ergodic com- 
ponents are maximal strongly connected subgraphs 
of the graph, but inverse is not true: there may ex- 
ist maximal strongly connected subgraphs that have 
outgoing edges and, therefore, are not sinks. 

We study ensembles of systems with a given 
graph and independent and well separated kinetic 
constants kij . This means that we study asymptotic 
behaviour of ensembles with independent identi- 
cally distributed constants, log-uniform distributed 
in sufficiently big interval log A; £ [a,/3], for a —^ 
—00, — > 00, or just a log- uniform distribution on 
infinite axis, log fc G M. 

4.1.2. Sinks and ergodicity 

If there is no other independent linear conserva- 
tion law, then the system is weakly ergodic. The 
weak ergodicity of the network follows from its topo- 
logical properties. 

The following properties are equivalent and each 
one of them can be used as an alternative definition 
of weak ergodicity: 

(i) There exist the only independent linear con- 
servation law for kinetic equations ([27]) (this 
is 6°(c) = Ci = const). 

(ii) For any normalized initial state c(0) {b^{c) = 
1) there exists a limit state 

c* = lim exp{Kt) c(0) 

t — >oo 

that is the same for all normalized initial con- 
ditions: For all c, 

lim exp(Xt)c = 6°(c)c*. 

t — *oo 

(iii) For each two vertices Ai, Aj [i ^ j) we can 
find such a vertex Ak that is reachable both by 
Ai and by Aj. This means that the following 
structure exists; 

Ai > . . . !■ Ak < . . . < Aj . 

One of the paths can be degenerated: it may 
he i — k or j — k. 

(iv) The network has only one minimal sink (one 
ergodic component). 
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For every monomolecular kinetic system, the Jor- 
dan cell for zero eigenvalue of matrix K is diagonal 
and the maximal number of independent linear con- 
servation laws (i.e. the geometric multiplicity of the 
zero eigenvalue of the matrix K) is equal to the max- 
imal number of disjoint ergodic components (mini- 
mal sinks). 

Let G = {Ai-^ , . . . } be an ergodic component. 
Then there exists a unique vector (normalized in- 
variant distribution) c*^ with the following proper- 
ties: cf = for i ^ {ii, . . . ii}, cf > for all i € 
{ii,...ii},b"{c^) = l,Kc(^ = 0. 

If Gi, . . . Gm are all ergodic components of the 
system, then there exist m independent positive lin- 
ear functionals b^{c), ... b"^{c) such that X]i=i ^* = 
6° and for each c 



lim exp{Kt)c 

t — >oo 



.Gi 



(30) 



So, for any solution of kinetic equations ((27|) . c{t), 
the limit at t — > oo is a linear combination of nor- 
malized invariant distributions c^' with coefficients 
&*(c(0)). In the simplest example, Ai ^ A2 ^ A3, 
Gi = {^i}, G2 = {A3}, components of vectors c"^^ , 
c'^^ are (1,0,0) and (0,0,1), correspondingly. For 
functionals h^'^ we get: 



&\c) =Cl + 



ki 



ki + k2 



C2; 6^(c) = ^C2+C3, 

ki + k2 

(31) 

where fci , k2 are rate constants for reaction A2 
Ai , and A2 ^ A3, correspondingly. We can mention 
that for well separated constants either fci ^ k2 or 
ki ^ /c2- Hence, one of the coefficients fci /(fci + fc2), 
fc2/(fci + fc2) is close to 0, another is close to 1. This 
is an example of the general zero-one law for multi- 
scale systems: for any I, i, the value of functional 6' 
((30)) on basis vector e*, 5'(e'), is either close to one 
or close to zero (with probability close to 1). 

We can understand better this asymptotics by us- 
ing the Markov chain language. For non-separated 
constants a particle in A2 has nonzero probability to 
reach Ai and nonzero probability to reach A3 . The 
zero-one law in this simplest case means that the dy- 
namics of the particle becomes deterministic: with 
probability one it chooses to go to one of vertices 
A2, A3 and to avoid another. Instead of branching, 
A2 Ai and A2 ^ A3, we select only one way: ei- 
ther ^2 ^ ^1 or ^ ^3- Graphs without branch- 
ing represent discrete dynamical systems. 



4.1.3. Decomposition of discrete dynamical systems 
Discrete dynamical system on a finite set V = 
{Ai,A2, . . . An} is a semigroup 1, (j), (j)^ , where (j) 
is a map (j) : V ^ V . Ai £ V is & periodic point, 
if (f)^{Ai) — Ai for some I > 0; else Ai is a tran- 
sient point. A cycle of period I is a sequence of I 
distinct periodic points A, 0(A), 0^(A), . . . (j)'-^^{A) 
with 0'(A) ^ A. A cycle of period one consists of 
one fixed point, (j){A) = A. Two cycles, C, C' either 
coincide or have empty intersection. 

The set of periodic points, V^, is always 
nonempty. It is a union of cycles: = UjCj. For 
each point A £ V there exist such a positive integer 
t{A) and a cycle C{A) = Cj that (/)«(A) e Q for 
q > t{A). In that case we say that A belongs to basin 
of attraction of cycle Cj and use notation Att{Cj) = 
{A I C{A) = Cj}. Of course, Q C Att{Cj). For dif- 
ferent cycles, Att{Cj)nAtt{Ci) = 0. If A is periodic 
point then t{A) = 0. For transient points t{A) > 0. 

So, the phase space V is divided onto subsets 
Att{Cj). Each of these subsets includes one cycle 
(or a fixed point, that is a cycle of length 1). Sets 
Att{Cj) are (/)-invariant: (/)(Att(Cj)) C Att{Cj). The 
set Att(Cj ) \ Cj consist of transient points and there 
exists such positive integer r that (/)'(Att(Cj)) = Cj 
ii q > T. 



4.2. Auxiliary discrete dynamical systems and 
relaxation analysis 

4.2.1. Auxiliary discrete dynamical system 

For each Ai, we define Kii clS the maximal kinetic 



constant for reactions Ai 



A, 



maxjlkji}. 



For correspondent j we use notation (t>(i): 4>{i) = 
argmax^{fcjj}. The function <j)[i) is defined under 
condition that for Ai outgoing reactions Ai Aj 
exist. Let us extend the definition: (f){i) = i if there 
exist no such outgoing reactions. 

The map (j) determines discrete dynamical system 
on a set of components V = {Ai}. We call it the 
auxiliary discrete dynamical system for a given net- 
work of monomolecular reactions. Let us decompose 
this system and find the cycles Cj and their basins 
of attraction, Att{Cj). 

Notice that for the graph that represents a dis- 
crete dynamic system, attractors are ergodic com- 
ponents, while basins are connected components. 

An auxiliary reaction network is associated with 
the auxiliary discrete dynamical system. This is the 



set of reactions Ai 



with kinetic constants 



Ki . The correspondent kinetic equation is 
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— —KiCi + ^ ^ 



(32) 



or in vector notations 



(33) 

For deriving of the auxiliary discrete dynamical 
system we do not need the values of rate constants. 
Only the ordering is important. Below we consider 
multiscale ensembles of kinetic systems with given 
ordering and with well separated kinetic constants 
(^^(i) ^ ^o-(2) ^ •■• for some permutation a). 

In the following, we analyze first the situation 
when the system is connected and has only one at- 
tractor. This can be a point or a cycle. Then, we 
discuss the general situation with any number of at- 
tractors. 



4.2.2. Eigenvectors for acyclic auxiliary kinetics 

Let us study kinetics for acyclic discrete dy- 
namical system (each vertex has one or zero outgo- 
ing reactions, and there are no cycles). Such acyclic 
reaction networks have many simple properties. For 
example, the nonzero eigenvalues are exactly minus 
reaction rate constants, and it is easy to find all left 
and right eigenvectors in explicit form. Let us find 
left and right eigenvectors of matrix K of auxiliary 
kinetic system (j32p for acyclic auxiliary dynamics. 
In this case, for any vertex there exists is an 
eigenvector. If Ai is a fixed point of the discrete dy- 
namical system (i.e. 4>{i) = i) then this eigenvalue is 
zero. If Ai is not a fixed point (i.e. 4>{i) ^ i and reac- 
tion Ai — > ^<t>{i) has nonzero rate constant k^) then 
this eigenvector corresponds to eigenvalue —Ki. For 
left and right eigenvectors of K that correspond to 
Ai we use notations P (vector-raw) and r* (vector- 
column), correspondingly, and apply normalization 
condition — II — 1. 

First, let us find the eigenvectors for zero eigen- 
value. Dimension of zero eigenspace is equal to the 
number of fixed points in the discrete dynamical sys- 
tem. If Ai is a fixed point then the correspondent 
eigenvalue is zero, and the right eigenvector has 
only one nonzero coordinate, concentration of Ai: 

To construct the correspondent left eigenvectors 
r for zero eigenvalue (for fixed point Ai), let us men- 
tion that could have nonzero value only if there 
exists such g > that — i (this q is unique 

because absence of cycles) . In that case (for q > 0), 



Hence, = ^^q), and ^ = 1 if ^^(j) = i for some 

For nonzero eigenvalues, right eigenvectors will be 
constructed by recurrence starting from the vertex 
Ai and moving in the direction of the flow. The con- 
struction is in opposite direction for left eigenvec- 
tors. Nonzero eigenvalues oi K (15^ are —K,i. 

For given i, is the minimal integer such that 
(j)^^{i) — 0'^*~'"^(z) (this is a relaxation time i.e. the 
number of steps to reach a fixed point). All indices 
{(f)^ {i) \ k = . . .Ti} are different. For right eigen- 
vector r* only coordinates r^^k^i^ {k = 0,1,... r^) 
could have nonzero values, and 

(is:r*)0fc+i(j) = -K0fc+i(j)r^fe+i(j) (i)^^'=(j) 
Hence, 

k 

r\l■^.^,,^ = r\k,,\ — 



^0H') i _ 1 r (') 



K^3 + l(i) 



3 = 
fc-1 

= — - — n 

(34) 

The last transformation is convenient for estimation 
of the product for well separated constants (compare 
to Q): 



+ - Ki I 0, if + < Kf, 

Kj _ J -1, if Ki > K^fe+l(j), 

K4,k + ll^^^ - Ki 1 0, if < 



(35) 



For left eigenvector P coordinate Z] could have 
nonzero value only if there exists such q > that 
0''(j) = * (this q is unique because the auxiliary 
dynamical system has no cycles). In that case (for 
'?>0), 



{l^K)j — —Kjlj + Kjl\,,s — —KiPj. 



Hence, 



q-l 

n K^kf^j) - Ki 



(36) 



For every fraction in (j36p the following estimate 
holds: 



K^k(^j) 



1, if K, 



> Ki 



0, if K0fc(j) < K 



(37) 
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As we can see, every coordinate of left and right 
eigenvectors of K (l34l) . ((36)) is either or ±1, or 
close to or to ±1 (with probability close to 1). We 
can write this asymptotic representation explicitly 
(analogously to dS])). For left eigenvectors, l\ = I 
and = 1 (for i ^ j) \i there exists such q that 
4>''{j) — i, and n^di^j-^ > Kj for all d = 0, . . . g — 
1, else Ij = 0. For right eigenvectors, rl — 1 and 
^^'=(i) = ~1 if '*0'=(i) < '^i and for all positive m < 
k inequality K0m(j) > Ki holds, i.e. k is first such 
positive integer that K^k < Ki (for fixed point Ap 
we use Kp = 0). Vector has not more than two 
nonzero coordinates. It is straightforward to check 
that in this asymptotic Pr^ = dij . 

In general, coordinates of eigenvectors Ij, are 
simultaneously nonzero only for one value j = i be- 
cause the auxiliary system is acyclic. On the other 
hand, /V^ = if i 7^ j, just because that are eigen- 
vectors for different eigenvalues, Ki and kj. Hence, 

For example, let us find the asymptotic of left and 
right eigenvectors for a branched acyclic system of 
reactions: 

Ai^A2^A3-^Ai^A5^As, A(i^A7^Ai 

where the upper index marks the order of rate con- 
stants: kq > K4 > Kr > K5 > K2 > K3 > Ki (Ki is 
the rate constant of reaction Ai — > ...). 

For zero eigenvalue, the left and right eigenvectors 
are 

= (1,1,1,1,1,1,1,1,1), ^ (0,0,0,0,0,0,0,1). 

For left eigenvectors, rows P, that correspond to 
nonzero eigenvalues we have the following asymp- 
totics: 

« (1,0,0,0,0,0,0,0), l"^ « (0,1,0,0,0,0,0,0), 
« (0,1, 1,0, 0,0, 0,0), (0,0,0,1,0,0,0,0), 
f « (0,0,0,1,1,1,1,0), « (0,0,0,0,0,1,0,0). 

f «:! (0,0,0,0,0,1,1,0) 

(38) 

For the correspondent right eigenvectors, columns 
r', we have the following asymptotics (we write 
vector-columns in rows): 

r^w (1,0, 0,0, 0,0, 0,-1), r2w(0, 1,-1, 0,0, 0,0,0), 
(0, 0, 1,0, 0,0, 0,-l),r4«(0, 0,0, 1,-1, 0,0,0), 
r5«(0,0,0,0,l,0,0,-l),r6«(0,0,0,0,0,l,-l,0), 
r7«(0, 0,0, 0,-1, 0,1,0). 

(39) 



4.2.3. The first case: auxiliary dynamical system is 
acyclic and has one attractor 

In the simplest case, the auxiliary discrete dynam- 
ical system for the reaction network W is acyclic and 
has only one attractor, a fixed point. Let this point 
be An {n is the number of vertices). The correspon- 
dent eigenvectors for zero eigenvalue are r" = Snj, 

— 1. For such a system, it is easy to find explicit 
analytic solution of kinetic equation p2[) . 

Acyclic auxiliary dynamical system with one at- 
tractor have a characteristic property among all aux- 
iliary dynamical systems: the stoichiometric vectors 
of reactions Ai — > ^^(i) form a basis in the subspace 
of concentration space with Ci = 0. Indeed, for 
such a system there exist n — 1 reactions, and their 
stoichiometric vectors are independent. On the other 
hand, existence of cycles implies linear connections 
between stoichiometric vectors, and existence of two 
attractors in acyclic system implies that the num- 
ber of reactions is less n — 1, and their stoichiometric 
vectors could not form a basis in n — 1-dimensional 
space. 

Let us assume that the auxiliary dynamical sys- 
tem is acyclic and has only one attractor, a fixed 
point. This means that stoichiometric vectors 7,/,(i)i 
form a basis in a subspace of concentration space 
with J2i Ci = 0. For every reaction Ai — > Ai the fol- 
lowing linear operators Qn can be defined: 



111 



'»/(70(p)p) = for j3 7^ i. (40) 



The kinetic equation for the whole reaction network 
(l28l) could be transformed in the form 



At 



1+ E 



3,1 {i^4>{j)) 



E 



hi (41) 



'-Qji Kc, 



where K is kinetic matrix of auxiliary kinetic equa- 
tion (|33p . By construction of auxiliary dynamical 
system, ku <^ Ki \i I ^ 4>{^)- Notice also that \Qji\ 
does not depend on rate constants. 

Let us represent system (|41l) in eigenbasis of K ob- 
tained in previous subsection. Any matrix B in this 
eigenbasis has the form B = {bij), bij — VBr^ = 
J2qs^q^gsri, where (bqs) is matrix B in the initial 
basis, l^ and r^ are left and right eigenvectors of K 



19 



In eigenbasis of K the Gershgorin esti- 
mates of eigenvalues and estimates of eigenvectors 
are much more efficient than in original coordinates: 
the system is stronger diagonally dominant. Trans- 
formation to this basis is an effective precondition- 
ing for perturbation theory that uses auxiliary ki- 
netics as a first approximation to the kinetics of the 
whole system. 

First of all, we can exclude the conservation law. 
Any solution of ((4T|) has the form c{t) ~ br" + c{t), 
where b = l"-c{t) = l"c{0) and J2i = 0. On the 
subspace of concentration space with J2i Q = we 
get 



dc 

— = (1 +£)diag{-Ki, 
at 



K„_l}c 



(42) 



where £ = (cij), \sij\ ^ 1, and diag{— — 
K„_i} is diagonal matrix with —ki, ... — on 
the main diagonal. If \eij \ <^ 1 then we can use the 
Gershgorin theorem and state that eigenvalues of 
matrix (1 -I- £)diag{— ki, ... — are real and 

have the form Ai — + Oi with \9i \ <ti Ki. 

To prove inequality \eij\ <C 1 (for ensembles with 
well separated constants, with probability close to 
1) we use that the left and right eigenvectors of K 
(I34p . (|36p are uniformly bounded under some non- 
degeneracy conditions and those conditions arc true 
for well separated constants. For ensembles with well 
separated constants, for any given positive g < 1 
and all i,j {i ^ j) the following inequality is true 



with probability close to 1: \Ki 



> QKi. Let us 



select a value of g and assume that this diagonal 
gap condition is always true. In this case, for every 
fraction in (l34l) , (l36t we have estimate 



1 

< -. 

9 



Therefore, for coordinates of right and left eigenvec- 
tors of if jSll), jSH) we get 



gq gn 



(43) 



We can estimate | and \Oi\/Ki from above as 
const X max/^0(5){fc;s/Ks}. So, the eigenvalues for 
kinetic matrix of the whole system (|4ip are real and 
close to eigenvalues of auxiliary kinetic matrix K 
(f33|l . For eigenvectors, the Gershgorin theorem gives 
no result, and additionally to diagonal dominance 
we must assume the diagonal gap condition. Based 
on this assumption, we proved the Gershgorin type 
estimate of eigenvectors in Appendix 1. In particu- 
lar, according to this estimate, eigenvectors for the 



whole reaction network are arbitrarily close to eigen- 
vectors of K (with probability close to 1). 

So, if the auxiliary discrete dynamical system is 
acyclic and has only one attractor (a fixed point), 
then the relaxation of the whole reaction network 
could be approximated by the auxiliary kinetics 

n-l 

cit) = (/"c(0))r" + ^(f c(0))r* exp(-K,t), (44) 

i=l 

For r and one can use exact formulas ([Ml) and 
(|36p or zero-one asymptotic representations based 
on (|37p and ([55]) for multiscale systems. This approx- 
imation (|44[) could be improved by iterative meth- 
ods, if necessary. 

4.2.4. The second case: auxiliary system has one 
cyclic attractor 

The second simple particular case on the way to 
general case is a reaction network with components 
Ai , . . . An whose auxiliary discrete dynamical sys- 
tem has one attractor, a cycle with period r > 1: 
An-T+i An-T+2 ■ ■ ■ An — > An-r+i (after 
some change of enumeration) . We assume that the 
limiting step in this cycle (reaction with minimal 
constant) is An — > An-r+i- If auxiliary discrete 
dynamical system has only one attractor then the 
whole network is weakly ergodic. But the attractor 
of the auxiliary system may not coincide with a sink 
of the reaction network. 

There are two possibilities: 

(i) In the whole network, all the outgoing reac- 
tions from the cycle have the form An-r+i 
An-T+j {i,j > 0). This means that the cycle 
vertices An-r+i, An-T+2, ■ ■ ■ An form a sink 
for the whole network. 

(ii) There exists a reaction from a cycle vertex 
An-T+i to Am, m < n — T. This means that 
the set {An-T+i, An-T+2, ■ ■ ■ An} is not a sink 
for the whole network. 

In the first case, the limit (for t — -> oo) distribution 
for the auxiliary kinetics is the well-studied station- 
ary distribution of the cycle An-r+i, An-T+2, ■ ■ ■ An 
described in Sec. MUM, dHl) (113, (fTSl). The set 
{An-T+i, An-T+2, ■ ■ ■ An} is the only ergodic com- 
ponent for the whole network too, and the limit dis- 
tribution for that system is nonzero on vertices only. 
The stationary distribution for the cycle An-r+i 
An-T+2 — > • . . An — > An-T+1 approximates the sta- 
tionary distribution for the whole system. To ap- 
proximate the relaxation process, let us delete the 
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limiting step An — > An-r+i from this cycle. By this 
deletion we produce an acyclic system with one fixed 
point, Am and auxiliary kinetic equation (|33p trans- 
forms into 

d ^ 

_ = Kqc = ^ tiCi70(i) (45) 

i=l 

As it is demonstrated, dynamics of this system ap- 
proximates relaxation of the whole network in sub- 
space ^^Ci = 0. Eigenvalues for are — (z < 
n) , the corresponded eigenvectors are represented by 
((34| . ([36)1 and zero-one multiscale asymptotic rep- 
resentation is based on (|57)l and (|35p . 
In the second case, the set 

An} 

is not a sink for the whole network. This means 
that there exist outgoing reactions from the cycle, 

An-r+i Aj with ^ {An-r+1, An—T+2, ■ ■ ■ An}- 

For every cycle vertex An-r+i the rate constant 
Kn~T+i that corresponds to the cycle reaction 
An-T+i ^ An-T+i+i IS much bigger than any other 
constant kj^ n-r+i that corresponds to a "side" reac- 
tion An-T+i ^ Aj (j n - T + i + 1): Kn-r+i > 

kj^ n-T+i- This is because definition of auxiliary dis- 
crete dynamical system and assumption of ensemble 
with well separated constants (multiscale asymp- 
totics) . This inequality allows us to separate motion 
and to use for computation of the rates of outgoing 
reaction An-r+i Aj the quasi steady state dis- 
tribution in the cycle. This means that we can glue 
the cycle into one vertex An_^^i with the corre- 
spondent concentration c\_^j^i — J2i<i<T Cn-r+i 
and substitute the reaction An-r+i ^ Aj by 
An_^_^i — !■ Aj with the rate constant renormal- 

ization: fcj,„_^+i = kJ^n-r+^c^lr+^/(^i-r+l■ By 
the superscript QS we mark here the quasistation- 
ary concentrations for given total cycle concentra- 
tion Cn_^^i. Another possibility is to recharge the 
link An-r+i ^ Aj to another vertex of the cycle 
(usually to An): we can substitute the reaction 
An-r+i Aj by the reaction An-r+q Aj with 
the rate constant renormalization: 

kj^n-r+q =^ kj_ n-r+iCn-r+i/ ^n-r+q- (46) 

The new rate constant is smaller than the initial one: 

kj^n-r+q < kj^n-r+i, bcCaUSe C^^^_,_^ < C^^r+q ^UC 

to definition. 

We apply this approach now and demonstrate 
its applicability in more details later in this sec- 
tion. For the quasi-stationary distribution on the 

cycle we get Cn-r+i = CnK.nl ^n-r+t (1 < « < t). 



The original reaction network is transformed by 
gluing the cycle {An-r+i, An-r+2, ■ ■ ■ An} into a 
point An_^_^_i. We say that components An-r+i, 
An-T+2, ■ ■ ■ An of the original system belong to 
the component An_^_^_l of the new system. All the 
reactions Ai — > Aj with i,j < n — t remain the 
same with rate constant kji . Reactions of the form 
Ai —^ Aj with i < n — T, j > n — t (incoming 
reactions of the cycle {An-r+i, An-r+2, ■ ■ ■ An}) 
transform into Ai — > Al^_^_^^ with the same rate 
constant kji. Reactions of the form Ai —> Aj with 
i > n — T, i < n — T (outgoing reactions of the 
cycle {An-r+i, An-r+2, ■ ■ ■ An}) transform into re- 
actions An_^^i — > Aj with the "quasistationary" 
rate constant k^^ = kjiKn/ Kn-r+i- After that, we 
select the maximal k'^f for given j: = 
maxi>„_,- fc^^. This kj^l^_^_^_^ is the rate constant 
for reaction An_^_^_i — s- Aj in the new system. Reac- 
tions Ai Aj with i,j > n T (internal reactions 
of the site) vanish. 

Among rate constants for reactions of the form 

An-r+i — > Am {m > TL ~ t) WC find 

K^i—r+i ^ "^&^{km,n — r+iKn / Kn-r+i} ■ (47) 
i.7n 

Let the correspondent «,to be ii, mi. 

After that, we create a new auxiliary discrete dy- 
namical system for the new reaction network on the 
set {Ai, . . . An-rT A\_^_^{}. We can describe this 
new auxiliary system as a result of transformation 
of the first auxiliary discrete dynamical system of 
initial reaction network. All the reactions from this 
first auxiliary system of the form Ai — s- Aj with 
*j j ^ n — T remain the same with rate constant Ki. 
Reactions of the form Ai —^ Aj with i < n -~ t, 
j > n — T transform into Ai —>■ A]^_^_^^ with the 
same rate constant k,;. One more reaction is to be 
added: An_^j^i A^i with rate constant Knlr+i- 
We "glued" the cycle into one vertex, and 
added new reaction from this vertex to A^i with 
maximal possible constant (|T7)) . Without this reac- 
tion the new auxiliary dynamical system has only 
one attractor, the fixed point An_^_^_i. With this ad- 
ditional reaction that point is not fixed, and a new 
cycle appears: Am^ ^ ... A}^-r+\ ^ • 

Again we should analyze, whether this new cycle 
is a sink in the new reaction network, etc. Finally, 
after a chain of transformations, we should come to 
an auxiliary discrete dynamical system with one at- 
tractor, a cycle, that is the sink of the transformed 
whole reaction network. After that, we can find sta- 
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Fig. 1. The main operation of the cycle surgery: on a step 
back we get a cycle Ai —>...—> Ar — » Ai with the limiting 
step At — » Ai and one outgoing reaction Ai Aj. We 
should delete the limiting step, reattach ("recharge") the 
outgoing reaction Ai — > Aj from Ai to At and change its 
rate constant k to the rate constant kkn^/ki. The new value 
of reaction rate constant is always smaller than the initial 
one: kkny^-^/ki < k if kn^ 7^ ki. For this operation only one 
condition A; <C fc; is necessary (k should be small with respect 
to reaction Ai A^+i rate constant, and can exceed any 
other reaction rate constant). 

tionary distribution by restoring of glued cycles in 
auxiliary kinetic system and applying formulas , 
(fT2ll (fT3|l . (fT5ll from Sec. [2 First, we find the station- 
ary state of the cycle constructed on the last itera- 
tion, after that for each vertex Aj that is a glued cy- 
cle we know its concentration (the sum of all concen- 
trations) and can find the stationary distribution, 
then if there remain some vertices that are glued cy- 
cles we find distribution of concentrations in these 
cycles, etc. At the end of this process we find all 
stationary concentrations with high accuracy, with 
probability close to one. 

As a simple example we use the following system, 
a chain supplemented by three reactions: 



Ae^A4, A5^A2, As^Ai, 



(48) 



where the upper index marks the order of rate con- 
stants. 

Auxiliary discrete dynamical system for the net- 
work (|48|) includes the chain and one reaction: 

A.l^A^^As^A^^A^^Ae^A^. 

It has one attractor, the cycle A4^A^-^Aq-^A4. 
This cycle is not a sink for the whole system, be- 
cause there exists an outgoing reaction A5^A2. 

By gluing the cycle A4^A5'^Aq-^A4 into a ver- 
tex A\ we get new network with a chain supple- 
mented by two reactions: 



^As^Al 



(49) 



Here the new rate constant is k. 



(k6 



(1) 

24 



k25Ke/n5 



k^Q is the limiting step of the cycle 
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Ai^A^^A^^Ai, K5 = h 

Here we can make a simple but important obser- 
vation: the new constant k\i^ — ^25^6/^5 has the 
same log-uniform distribution on the whole axis as 
constants /C25, and K5 have. The new constant k^^ 
depends on ^25 and the internal cycle constants Kg 
and /C5, and is independent from other constants. 
Of course, k24 < K5, but relations between ^24'' 
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and fci3 are a priori unknown. Both ordcrings, k'^^ 
fci3 and fc24'' < A;i3, are possible, and should be con- 
sidered separately, if necessary. But for both ordcr- 
ings the auxiliary dynamical system for network (|49p 
is 

A^^A2^A^^A\^A2 

(of course, k^^' < K3 < K2 < ki). It has one at- 
tractor, the cycle A2^A^^A\^A2. This cycle is 
not a sink for the whole system, because there ex- 
ists an outgoing reaction A^~^Ai. The limiting con- 
stant for this cycle is ^4"'^'' — k^24 — ^25^46/^65- We 
glue this cycle into one point, A\. The new trans- 
formed system is very simple, it is just a two step 
cycle: Ai^A^^Ai. The new reaction constant is 

^12'' = fci3'«4^V'^3 = ^13^5^46/(^65^43)- The aux- 
iliary discrete dynamical system is the same graph 
Ai^A\^Ai , this is a cycle, and we do not need fur- 
ther transformations. 

Let us find the steady state on the way back, from 
this final auxiliary system to the original one. For 
steady state of each cycle we use formula ([TB]) . 

The steady state for the final system is ci = 



bk^^2 1^21, C2 = 6(1 — k\^2 The component A2 

includes the cycle A2^Aj,^A\^A2- The steady 
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state of this cycle is C2 



c^'^h^^^ Ik., r^^^ 



(1) 



2 ^24 /'^32, C3 



32 ^ 



.247^43). The 



component Al includes the cycle A4^Ar,-^Ae-^A4. 



The steady state of this cycle is C4 



^46/^54, 

C5 = ^fc46/fc65, C6 = 4^''(1 - k4e/k54 - /c46/^65)- 

For one catalytic cycle, relaxation in the subspace 
Ci = is approximated by relaxation of a chain 
that is produced from the cycle by cutting the lim- 
iting step (Sec.[2]). For reaction networks under con- 
sideration (with one cyclic attractor in auxiliary dis- 
crete dynamical system) the direct generalization 
works: for approximation of relaxation in the sub- 
space J2i Ci = it is sufficient to perform the fol- 
lowing procedures: 

~ To glue iteratively attractors (cycles) of the aux- 
iliary system that are not sinks of the whole sys- 
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tern; 

- To restore these cycles from the end of the first 
procedure to its beginning. For each of cycles (in- 
cluding the last one that is a sink) the limited 
step should be deleted, and the outgoing reaction 
should be reattached to the head of the limiting 
steps (with the proper normalization), if it was 
not deleted before as a limiting step of one of the 
cycles. 

The heads of outgoing reactions of that cycles 
should be reattached to the heads of the limiting 
steps. Let for a cycle this limiting step be Am Aq. 
If for a glued cycle A*' there exists an outgoing reac- 
tion A^ Aj with the constant k (|47)) . then after 
restoration we add the outgoing reaction Am — > Aj 
with the rate constant k. Kinetic of the resulting 
acyclic system approximates relaxation of the initial 
network (under assumption of well separated con- 
stants, for given ordering, with probability close to 

Let us construct this acyclic network for the same 
example dH]). The final cycle is Ai^Al^Ai. The 
limiting step in this cycle is A^-^Ai. After cut- 
ting we get Ai^A\. The component A\ is glued 
cycle A2^A'i^A\^A2. The reaction Ai^Al cor- 
responds to the reaction Ai^A2 (in this case, this 
is the only reaction from Ai to cycle; in other 
case one should take the reaction from Ai to cycle 
with maximal constant). The limiting step in the 
cycle is A\^A2. After cutting, we get a system 
Ai^A2^A3^Al. The component Al is the glued 
cycle A^^A^^Aq^A^ from the previous step. The 
limiting step in this cycle is Aq-^A4. After restoring 
this cycle and cutting the limiting step, we get an 
acyclic system Ai^A2-^A3^A4^A5^Aq (as one 
can guess from the beginning: this coincidence is 
provided by the simple constant ordering selected 
in (j48|) ). Relaxation of this system approximates 
relaxation of the whole initial network. 

To demonstrate possible branching of described 
algorithm for cycles surgery (gluing, restoring and 
cutting) with necessity of additional orderings, let 
us consider the following system: 

Ai^A2^A3^A4^A5^A3, A4^A2, (50) 

The auxiliary discrete dynamical system for reac- 
tion network ([50| is 

A,1.A2^A3^A4^A,^A3. 

It has only one attractor, a cycle A3^A4-^A5^A3. 
This cycle is not a sink for the whole network (l50|) be- 



cause reaction A4-^A2 leads from that cycle. After 
gluing the cycle into a vertex A\ we get the new net- 
work Ai^A2^Al^A2. The rate constant for the 
reaction ^3— S-A2 is k^^ = A:24^35/fc54, where kij is 
the rate constant for the reaction Aj Ai in the 
initial network (/cas is the cycle limiting reaction). 
The new network coincides with its auxiliary system 
and has one cycle, A2^Al^A2. This cycle is a sink, 
hence, we can start the back process of cycles restor- 
ing and cutting. One question arises immediately: 
which constant is smaller, k32 or fcjs- The smallest 
of them is the limiting constant, and the answer de- 
pends on this choice. Let us consider two possibili- 
ties separately: (1) fc32 > and (2) ^32 < fc23. Of 
course, for any choice the stationary concentration 
of the source component Ai vanishes: ci = 0. 

(1) Let as assume that A:32 > fc23- this case, 
the steady state of the cycle A2^Al^A2 is (ac- 
cording to (US])) C2 = ^'^23/^32, C3 = 6(1 - fc23/fc32), 

where b = ^Ci. The component A^ is a glued 
cycle A3^A4-^A5^A3. Its steady state is C3 = 
4^35/^43, C4 = clkss/k^i, C5 = €3(1 - A;35/fc43 - 
^35/^54)- 

Let us construct an acyclic system that approx- 
imates relaxation of (|50p under the same assump- 
tion (1) /c32 > ^23. The final auxiliary system after 
gluing cycles is Ai^A2-^Al-^A2. Let us delete the 
limiting reaction ^13^^12 from the cycle. We get an 
acyclic system Ai^A2^A^. The component ^3 is 
the glued cycle A3-^A4^A5^A3. Let us restore this 
cycle and delete the limiting reaction A^^A^. We 
get an acyclic system Ai^^2-^^3^^4^^5- Re- 
laxation of this system approximates relaxation of 
the initial network (|50p under additional condition 
/C32 > fc^3. 

(2) Let as assume now that /c32 < ^23- this case, 
the steady state of the cycle ^12-^^3 -^^2 is (accord- 
ing to C2 = b{l - k32/kl3), cl = bk32/kls. The 
further analysis is the same as it was above: C3 = 
4^35/^43, C4 = 4^:35/^54, C5 = 4(1 - A;35/fc43 - 
^35/^54) (with another C3). 

Let us construct an acyclic system that approxi- 
mates relaxation of (f50|) under assumption (2) ^32 < 
/C23 ■ The final auxiliary system after gluing cycles is 
the same, Ai^A2-^Al^A2, but the limiting step 
in the cycle is different, ^12-^^3. After cutting this 
step, we get acyclic system Ai^A2^Al, where the 
last reaction has rate constant ^23 ■ 

The component ^3 is the glued cycle 

A3^A4^A5^A3 . 

Let us restore this cycle and delete the limiting re- 
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action A^^A^. The connection from glued cycle 
A^^A2 with constant transforms into connec- 
tion A5^A2 with the same constant ^23. 
We get the acyclic system: 

The order of constants is now known: k2i > k^s > 
> fc23, and we can substitute the sign "?" by 

"4":^3^A4^A5^A2. 

For both cases, fc32 > k^^ (^23 — ^24^35/^54) and 
^32 < ^23 it is easy to find the eigenvectors explicitly 
and to write the solution to the kinetic equations in 
explicit form. 

4.3. The general case: cycles surgery for auxiliary 
discrete dynamical system with arbitrary family of 
attractors 

In this subsection, we summarize results of relax- 
ation analysis and describe the algorithm of approx- 
imation of steady state and relaxation process for 
arbitrary reaction network with well separated con- 
stants. 

4.3.1. Hierarchy of cycles gluing 

Let us consider a reaction network W with a given 
structure and fixed ordering of constants. The set of 
vertices of W is and the set of elementary reactions 
is TZ. Each reaction from TZ has the form Ai Aj, 
Ai,Aj E A. The correspondent constant is kji. For 
each Ai E A we define k,; — maxjjfcji} and = 
argmaXj{fcji}. In addition, = i if kji = for all 

The auxiliary discrete dynamical system for the 
reaction network W is the dynamical system $ = 
$w defined by the map on the set A. Auxiliary 
reaction network V = Vyy has the same set of ver- 
tices A and the set of reactions Ai — > ^^(i) with re- 
action constants Ki. Auxiliary kinetics is described 
by c = Kc, where Kij = —KjSij + KjSi^(^j-j. 

Every fixed point of $vv is also a sink for the reac- 
tion network W. If all attractors of the system <I>w 
are fixed points Afi,Af2,--.EA then the set of sta- 
tionary distributions for the initial kinetics as well 
as for the auxiliary kinetics is the set of distributions 
concentrated the set of fixed points {Afi,Af2, ■■■}■ 
In this case, the auxiliary reaction network is acyclic, 
and the auxiliary kinetics approximates relaxation 
of the whole network W. 

In general case, let the system $vv have sev- 
eral attractors that are not fixed points, but cycles 



Ci,C2,... with periods ri,T2,... > 1. By gluing 
these cycles in points, we transform the reaction 
network W into W^. The dynamical system ^vv 
is transformed into For these new system and 
network, the connection <I>^ = <I>yyi persists: is 
the auxiliary discrete dynamical system for W^. 

For each cycle, Q, we introduce a new vertex A\ 
The new set of vertices, A^ = AU {A^^A^, ...} \ 
{UiCi) (we delete cycles Ci and add vertices ^'). 

All the reaction between A B [A, B E A) can 
be separated into 5 groups: 

(i) both A,B ^ U-jCi; 
(n) A ^ UiCi, but B e Ci] 
(in) A e Ci, but B ^ U^Q; 

(iv) A&C,,BEC,,i^r, 

(v) A,BE C,. 

Reactions from the first group do not change. Reac- 
tion from the second group transforms into A A'' 
(to the whole glued cycle) with the same constant. 
Reaction of the third type changes into A^ ^ B with 
the rate constant renormalization ([46]): let the cy- 
cle C" be the following sequence of reactions Ai 
A2 — !■ ---Ar- Ai, and the reaction rate constant 
for Ai Ai+i is ki {k^ for At^ Ai). For the lim- 
iting reaction of the cycle C'i we use notation fciim i. 
If A = Aj and k is the rate reaction for A ^ B, 
then the new reaction A'^ ^ B has the rate con- 
stant /c/ciini i/kj. This corresponds to a quasistation- 
ary distribution on the cycle ^3]) . It is obvious that 
the new rate constant is smaller than the initial one: 
^^lim i/kj < k, because /ciim i < kj due to definition 
of limiting constant. The same constant renormal- 
ization is necessary for reactions of the fourth type. 
These reactions transform into A* A^ . Finally, 
reactions of the fifth type vanish. 

After we glue all the cycles of auxiliary dynami- 
cal system in the reaction network W, we get W^. 
Strictly speaking, the whole network is not nec- 
essary, and in efficient realization of the algorithm 
for large networks the computation could be signifi- 
cantly reduced. What we need, is the correspondent 
auxiliary dynamical system $^ = <&vvi with auxil- 
iary kinetics. 

To find the auxiliary kinetic system, we should 
glue all cycles in the first auxiliary system, and then 
add several reactions: for each A* it is necessary to 
find in the reaction of the form A^ ^ B with 
maximal constant and add this reaction to the aux- 
iliary network. If there is no reaction of the form 
A* — !■ S for given i then the point A^ is the fixed 
point for and vertices of the cycle Ci form a sink 
for the initial network. 



24 



After that, we decompose the new auxihary dy- 
namical system, find cycles and repeat gluing. Ter- 
minate when all attractors of the auxiliary dynami- 
cal system $™ become fixed points. 

4.3.2. Reconstruction of steady states 

After this termination, we can find all steady 
state distributions by restoring cycles in the auxil- 
iary reaction network V™. Let ^/2i •■•be fixed 
points of The set of steady states for V™ is 
the set of all distributions on the set of fixed points 
{AJ\, AJ\, ...}. Let us take one of these distribu- 
tions, c = {cj\, c™2, •■•) (we mark the concentrations 
by the same indexes as the vertex has; other Ci = 0). 

To make a step of cycle restoration we select those 
vertexes AJl that are glued cycles and substitute 
them in the list A'J\,Aj2, ■■■ by all the vertices of 
these cycles. For each of those cycles we find the lim- 
iting rate constant and redistribute the concentra- 
tion between the vertices of the correspondent 
cycle by the rule (I13|) (with & = c^^). As a result, we 
get a set of vertices and a distribution on this set of 
vertices. If among these vertices there are glued cy- 
cles, then we repeat the procedure of cycle restora- 
tion. Terminate when there is no glued cycles in the 
support of the distribution. The resulting distribu- 
tion is the approximation to a steady state of W, 
and all steady states for W can be approximated by 
this method. 

In order to construct the approximation to the 
basis of stationary distributions of W, it is sufficient 
to apply the described algorithm to distributions 
concentrated on a single fixed point AJ-^, Cj" = Sij, 
for every j. 

The steady state approximation on the base of the 
rule (|13p is a linear function of the restored-and-cut 
cycles rate limiting constants. It is the first order 
approximation. 

The zero order approximation also makes sense. 
For one cycle is gives (|14|) : all the concentration is 
collected at the start of the limiting step. The algo- 
rithm for the zero order approximation is even sim- 
pler than for the first order. Let us start from the dis- 
tributions concentrated on a single fixed point AJ^^, 
c^j — Sij for some i. If this point is a glued cycle 
then restore that cycle, and find the limiting step. 
The new distribution is concentrated at the starting 
vertex of that step. If this vertex is a glued cycle, 
then repeat. If it is not then terminate. As a result 
we get a distribution concentrated in one vertex of 
A. 



4.3.3. Dominant kinetic system for approximation 
of relaxation 

To construct an approximation to the relaxation 
process in the reaction network W, we also need to 
restore cycles, but for this purpose we should start 
from the whole glued network V™ on A™ (not only 
from fixed points as we did for the steady state ap- 
proximation). On a step back, from the set to 
A™'~^ and so on some of glued cycles should be re- 
stored and cut. On each step we build an acyclic re- 
action network, the final network is defined on the 
initial vertex set and approximates relaxation of W. 

To make one step back from V™ let us select the 
vertices of A™ that are glued cycles from V™~^. Let 
these vertices be A™, A™, .... Each ^™ corresponds 
to a glued cycle from V"-i, A™"^ A^"^ ^ 
...A""^ ^ of the length n. We assume that 

the limiting steps in these cycles are A^~^ ^Ii"^ ■ 
Let us substitute each vertex A™ in V™ by ver- 
tices A™-\ A;^-!, ...A™-i and add to V" reactions 
A""^ ^ A™-i ^ -A^r^ (tliat are the cycle reac- 
tions without the limiting step) with correspondent 
constants from V™~^. 

If there exists an outgoing reaction A™ B in 
then we substitute it by the reaction A™^^ 
B with the same constant, i.e. outgoing reactions 
A,™ ... are reattached to the heads of the limiting 
steps. Let us rearrange reactions from V" of the 
form B A™. These reactions have prototypes in 
ym-i (^j-jgfQj-g tjig lagi; gluing). We simply restore 

these reactions. If there exists a reaction A™ A™ 
then we find the prototype in V""~^, A ^ B, and 
substitute the reaction by A™^^ — > B with the same 
constant, as for A™ A™. 

After that step is performed, the vertices set is 
^ but the reaction set differs from the reactions 
of the network V™"^: the limiting steps of cycles are 
excluded and the outgoing reactions of glued cycles 
are included (reattached to the heads of the limiting 
steps). To make the next step, we select vertices of 
A™^^ that are glued cycles from V""~^, substitute 
these vertices by vertices of cycles, delete the lim- 
iting steps, attach outgoing reactions to the heads 
of the limiting steps, and for incoming reactions re- 
store their prototypes from V"*"^, and so on. 

After all, we restore all the glued cycles, and con- 
struct an acyclic reaction network on the set A. This 
acyclic network approximates relaxation of the net- 
work W. We call this system the dominant system 
of W and use notation dommod(yV). 
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Fig. 2. Gluing of cycles for the prism of reactions with 
a given ordering of rate constants in the case of two at- 
tractors in the auxiliary dynamical system: (a) initial re- 
action network, (b) auxiliary dynamical system that con- 
sists of two cycles, (c) connection between cycles, (d) glu- 
ing cycles into new components, (e) network with 
glued vertices, (f) an example of dominant system in 
^41^32/^21 and > fc^fj (by 

max{fc4ifc32/A:2i, H2, fc63*:32/fcl3} and 
k^ek^A/k/^o), the order of constants in the dominant 
system is: k2\ > kif, > fcia > fces > ^41^32/^21- 



the case when fcj^ 
definition, k^^ 

"12 



4.4. Example: a prism of reactions 

Let us demonstrate work of the algorithm on a 
typical example, a prism of reaction that consists of 
two connected cycles fFig. l2l3p . Such systems appear 
in many areas of biophys ics and biochemis try (see, 
for example, the paper of lKurzvnski 

For the first example we use the reaction rate con- 
stants ordering presented in Fig. [5^. For this order- 
ing, the auxiliary dynamical system consists of two 
cycles (Fig. [2Jd) with the limiting constants fc54 and 
^32 , correspondingly. These cycles are connected by 
four reaction (Fig. [21;). We glue the cycles into new 
components A} and A2 (Fig. and the reaction 
network is transformed into A\ <~* A^. Following the 
general rule (fc^ — fc/cum/fej), we determine the rate 
constants: for reaction A\ A2 

= max{fc4ifc32/fc2i, fc52, fc63fc32/fci3}, 

A\ 

There are six possible orderings of the constant 
combinations: three possibilities for the choice of 
and for each such a choice there exist two possibili- 

trlGS . ^12 ^21 ^ ^12 " 

The zero order approximation of the steady state 
depends only on the sign of inequality between k\^ 



and for reaction A\ 



^12 = fc.'H 



and k\2. If k\i 3> fci2 then almost all concentration 
in the steady state is accumulated inside A\. After 
restoring the cycle A/^ ^ A^ ^ Aq ^ A4 we find 
that in the steady state almost all concentration is 
accumulated in A4 (the component at the beginning 
of the limiting step of this cycle, A4 A^). Finally, 
the eigenvector for zero eigenvalue is estimated as 
the vector column with coordinates (0, 0, 0, 1, 0, 0). 

If, inverse, k^-^ <C fcjj then almost all concentra- 
tion in the steady state is accumulated inside Al . Af- 
ter restoring the cycle Ai ^ A2 ^ A^ ^ Ai we find 
that in the steady state almost all concentration is 
accumulated in A2 (the component at the beginning 
of the limiting step of this cycle, A2 A3). Finally, 
the eigenvector for zero eigenvalue is estimated as 
the vector column with coordinates (0, 1, 0, 0, 0, 0). 

Let us find the first order (in rate limiting con- 
stants) approximation to the steady states. If fc^i ^ 
^22 then kl2 is the rate limiting constant for the cycle 
A\ A2 and almost all concentration in the steady 
state is accumulated inside A\: c\ ~ 1 — fci2/^2i 
and c\ K, k\2/k\-^. Let us restore the glued cycles 
(Fig. (2). In the upper cycle the rate limiting con- 
stant is fc32, hence, in steady state almost all concen- 
tration of the upper cycle, c}, is accumulated in A2: 
C2 « c\{l - k32/ki3 - ^32/^21), C3 « c\k32/ki3, and 
ci « c\k32/k2i- In the bottom cycle the rate limiting 
constant is ^54, hence, C4 « 02(1—^54/^65 — ^54/^:46), 
C5 « clk54/ke5 and cq ~ clk^i/k^e. 

If, inverse, k^i <SC kl2 then k^i is the rate limiting 
constant for the cycle Al ^ A\ and almost all con- 
centration in the steady state is accumulated inside 
A}: c} ~ 1 — ^21/^12 and c\ ~ k\i/k\2. For distribu- 
tions of concentrations in the upper and lover cycles 
only the prefactors c}, c\ change their values. 

For analysis of relaxation, let us analyze one of 
the six particular cases separately. 

1. k\i = ^41^32/^21 and fcji > ^12- 

In this case, the finite acyclic auxiliary dynamical 
system, <i)™ = $1, is A\ 
constant k],^ = ^41/032/^21, and is A\ 
restore both cycles and delete the limiting reactions 
A2 A3 and A4 A5. This is the common step 
for all cases. Following the general procedure, we 
substitute the reaction Al A\ by A2 A4 with 
the rate constant ^21 = ^41^32/^21 (because A2 is 
the head of the limiting step for the cycle Ai 
A2 A3 ^ Al, and the prototype of the reaction 
Al A\ is in that case Ai ^ A4. 

We find the dominant system for relaxation de- 
scription: reactions A3 Ai A2 and A^ 
Aq A4 with original constants, and reaction 



A2 with reaction rate 
Al We 
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A2 A4 with the rate constant fcjj^ = A:4i/c32/fc2i. 

This dominant system graph is acycUc and, more- 
over, represents a discrete dynamical system, as it 
should be (not more than one outgoing reaction for 
any component). Therefore, we can estimate the 
eigenvalues and eigenvectors on the base of formulas 
([55)) . ([57)1 . It is easy to determine the order of con- 
stants because = ^41^32 7^21: this constant is the 
smallest nonzero constant in the obtained acyclic 
system. Finally, we have the following ordering of 
constants: A3^Ai^A2'^A4, A^^Ae^A^. 

So, the eigenvalues of the prism of reaction for 
the given ordering are (with high accuracy, with 
probability close to one) — fc2i < —k^e < — fcia < 
— ^65 < —^41^32/^21- The relaxation time is r « 
^21/(^41^32)- 

We use the same notations as in previous sections: 
eigenvectors and correspond to the eigenvalue 
— Ki, where Ki is the reaction rate constant for the 
reaction Ai .... The left eigenvectors are: 

(1,0,0,0,0,0), (1,1,1,0,0,0), 

w (0,0,1,0,0,0), l"^ w (1,1,1,1,1,1), (51) 

« (0,0,0,0,1,0), « (0,0,0,0,0,1). 

The right eigenvectors r* are (we represent vector 
columns as rows): 

(1,-1,0,0,0,0), (0,1,0,-1,0,0), 

w (0, -1, 1, 0, 0, 0), r"^ w (0, 0, 0, 1, 0, 0), (52) 

(0,0,0,-1,1,0), rf^« (0,0,0,-1,0,1) 

The vertex A4 is the fixed point for the discrete dy- 
namical system. There is no reaction A4 .... For 
convenience, we include the eigenvectors and r** 
for zero eigenvalue, K4 = 0. These vectors corre- 
spond to the steady state: is the steady state vec- 
tor, and the functional is the conservation law. 

The correspondent approximation to the general 
solution of the kinetic equation for the prism of re- 
action (Fig. Oi) is: 

6 

c{t) = ^'(''' c(0)) exp(-K,t). (53) 

i=l 

Analysis of other five particular cases is similar. Of 
course, some of the eigenvectors and eigenvalues can 
differ. 

Of course, different ordering can lead to very dif- 
ferent approximations. For example, let us consider 
the same prism of reactions, but with the ordering 
of constants presented in Fig. [3^. The auxiliary dy- 
namical system has one cycle (Fig.[3)D) with the lim- 
iting constant fcsg. This cycle is not a sink to the 
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Fig. 3. Gluing of a cycle for the prism of reactions with a given 
ordering of rate constants in the case of one attractors in the 
auxiliary dynamical system: (a) initial reaction network, (b) 
auxiliary dynamical system that has one attractor, (c) outgo- 
ing reactions from a cycle, (d) gluing of a cycle into new com- 
ponent, (e) network with glued vertices, (f) an example 
of dominant system in the case when = kiQ, and, there- 
fore > ksi (by definition, k^ = max{fc4ifc36/fc2li '^46}); 
this dominant system is a linear chain that consists of some 
reactions from the initial system (no nontrivial monomials 
among constants). Only one reaction rate constant has in 
the dominant system new number (number 5 instead of 9). 

initial network, there are outgoing reactions from 
its vertices (Fig. [5};) . After gluing, this cycles trans- 
forms into a vertex A\ (Fig. [3Ji). The glued net- 
work, (Fig- [3fe), lias two vertices, A4 and A\ 
the rate constant for the reaction A4 A\ is fc54, 
and the rate constant for the reaction A\ A4 
is — max{fc4iA:36/fc2i, fc4g}. Hence, there are not 
more than four possible versions: two possibilities 
for the choice of k^ and for each such a choice there 
exist two possibilities: k^ > A:54 or k^ < ^54 (one of 
these four possibilities cannot be realized, because 
k46 > ^54))- 

Exactly as it was in the previous example, the 
zero order approximation of the steady state de- 
pends only on the sign of inequality between k^ and 
k54. If k^ <C fc54 then almost all concentration in the 
steady state is accumulated inside A^. After restor- 
ing the cycle A3 ^ Ai ^ A2 ^ A^ ^ Aq ^ A3 we 
find that in the steady state almost all concentration 
is accumulated in Aq (the component at the begin- 
ning of the limiting step of this cycle, Aq — > ^3). 
The eigenvector for zero eigenvalue is estimated as 
the vector column with coordinates (0, 0, 0, 0, 0, 1). 

If k^ 3> k54 then almost all concentration in the 
steady state is accumulated inside A'^ . This vertex is 
not a glued cycle, and immediately we find the ap- 
proximate eigenvector for zero eigenvalue, the vec- 
tor column with coordinates (0, 0, 0, 1, 0, 0). 
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Let us find the first order (in rate limiting con- 
stants) approximation to the steady states. If <C 
^54 then is the rate hmiting constant for the cycle 
A\ ^ A4 and almost all concentration in the steady 
state is accumulated inside A}: « 1 — fc^/fc54 and 
C4 « k^/k54. Let us restore the glued cycle (Fig.[3|). 
The limiting constant for that cycle is fcse, cq « 
c}(l - kse/kis - fcae/fei - kz&/k^2 - fcse/fces), C3 « 
c}fc36/fci3, ci Pi c\kz&/k2i, C2 « c\k-iQ/k^2, and C5 « 

If k^ 3> A;54 then fc54 is the rate limiting constant 
for the cycle A\ ^ A4 and almost all concentra- 
tion in the steady state is accumulated inside A^: 
C4 ~ 1 — k^i/k^ and c\ ~ kc^^/k^ . In distribution of 
concentration inside the cycle only the prefactor c\ 
changes. 

Let us analyze the relaxation process for one of the 
possibilities: fc^ = k^Q, and, therefore k^ > k^4. We 
restore the cycle, delete the limiting step, transform 
the reaction A\ —> A4 into reaction Aq —> A4 with 
the same constant k^ = ^45 and get the chain with 
ordered constants: A3^Ai^A2^A5^Ae-^A4. 
Here the nonzero rate constants kij have the same 
value as for the initial system (Fig. [3^). The re- 
laxation time is T w 1/^:46. Left eigenvectors are 
(including for the zero eigenvalue): 

f « (1,0,0,0,0,0), P « (1,1,1,0,0,0), 

(0,0,1,0,0,0); (1,1,1,1,1,1), (54) 
(0,0,0,0,1,0), (1,1,1,0,1,1). 

Right eigenvectors are (including for the zero 
eigenvalue) : 

w (1,-1,0,0,0,0), « (0,1,0,0,0,-1), 
r3 « (0, ~1, 1, 0, 0, 0), « (0, 0, 0, 1, 0, 0), (55) 
(0,0,0,0,1,-1), rf^« (0,0,0,-1,0,1). 

Here we represent vector columns as rows. 

For the approximation of relaxation in that order 
we can use (1551) . 

5. The reversible triangle of reactions: the 
simple example case study 

In this section, we illustrate the analysis of dom- 
inant systems on a simple example, the reversible 
triangle of reactions. 



Ai ^ A2 ^ A-s ^ Ai . 



(56) 




Fig. 4. Four possible auxiliary dynamical systems for the re- 
versible triangle of reactions with k2i > kij for ^ (2, 1): 
(a) ki2 > kz2, k23 > fcis; (b) ki2 > fe32, ki3 > k23; (c) 
^32 > fci2, ^23 > fcl3; (d) fc32 > ^12, ^13 > *:23- For each 
vertex the outgoing reaction with the largest rate constant 
is represented by the solid bold arrow, and other reactions 
are represented by the dashed arrows. The digraphs formed 
by solid bold arrows are the auxiliary discrete dynamical 
systems. Attractors of these systems are isolated in frames. 

triangle ([SS]) is not obligatory a closed system. We 
can assume that it is a subsystem of a larger sys- 
tem, and any reaction Ai — > Aj represents a reac- 



tion of the form , 



-Ai A, 



, where unknown 



This triangle appeared in many works as an ideal 
object f or a case study. Our favorite example is the 
work of IWei fc Prateii (jl962l ). Now in our study the 



but slow components are substituted by dots. This 
means that there are no obligatory relations between 
reaction rate constants, first of all, no detailed bal- 
ance relations, and six reaction rate constants are 
arbitrary nonnegative numbers. 

There exist 6! = 720 orderings of six reaction rate 
constants for this triangle, but, of course, it is not 
necessary to consider all these orderings. First of 
all, because of the permutation symmetry, we can 
select an arbitrary reaction as the fastest one. Let 
the reaction rate constant ^21 for the reaction Ai — > 
A2 is the largest. (If it is not, we just have to change 
the enumeration of reagents.) 

First of all, let us describe all possible auxiliary dy- 
namical systems for the triangle (|56p . For each ver- 
tex, we have to select the fastest outgoing reaction. 
For Ai, it is always Ai ^ A2, because of our choice 
of enumeration (the higher scheme in Fig. . There 
exist two choices of the fastest outgoing reaction for 
two other vertices and, therefore, only four versions 
of auxiliary dynamical systems for (I56|) (Fig. [4|) . 

Because of the choice of enumeration, the vec- 
tors of logarithms of reaction rate constants form 
a convex cone in which is described by the sys- 
tem of inequalities lnA:2i > Infcy, 7^ (2,1)- 
For each of the possible auxiliary systems (Fig |4]) 
additional inequalities between constants should 
be valid, and we get four correspondent cones in 

. This cones form a partitions of the initial one 
(we neglect intersections of faces which have zero 
measure). Let us discuss the typical behaviour of 
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systems from this cones separately. (Let us remind 
that if in a cone for some values of coefficients 6ij , 
Cii JZij^^ij^^f'ij < SijCijln^, then, typically in 
this cone J^ij % < K + J2xj Qj In hj for any 
positive K. This means that typically J| k^-^ <^ 

Uij ^fj' ■) 

5.1. Auxiliary system (a): Ai <-> <— A3; 

kl2 > k32, k23 > ki3 

5.1.1. Gluing cycles 

The attractor is a cycle (with only two vertices) 
<-!■ A2. This is not a sink, because two outgoing 
reactions exist: Ai A3 and A2 — > A3. They are 
relatively slow: fcai <C k2i and fc32 ^ /ci2 • The limit- 
ing step in this cycle is A2 Ai with the rate con- 
stant ki2. We have to glue the cycle Ai ^ A2 into 
one new component Al and to add a new reaction 
A\ — > A3 with the rate constant 



fcgi = max{fc32, k3iki2/k2i} ■ 



(57) 



This is a particular case of (l46|) . (f47|) . 

As a result, we get a new system, A\ ^ A3 with 
reaction rate constants (for A J A3) and initial 
^23 (for Al ^ A3). This cycle is a sink, because 
it has no outgoing reactions (the whole system is a 
trivial example of a sink) . 

5.1.2. Steady states 

To find the steady state, we have to compute the 
stationary concentrations for the cycle A\ ^ A3, c\ 
and C3. We use the standard normalization condition 
c} -I- C3 — 1 . On the base of the general formula for 
a simple cycle pT|) we obtain: 



1 



1 

fc23 



'''31 



-,^3 = ^. (58) 



After that, we can calculate the concentrations of 
Al and A2 with normalization C1+C2 = c\. Formula 
([TT|) gives: 



1 

'£21 



— > C2 = — . (59) 
ki2 



We can simplify the answer using inequalities be- 
tween constants, as it was done in formulas (|12p . 
(fTS]) . For example, + 17^ ~ 1771 because fei ^ 
ki2. It is necessary to stress that we have used the 
inequalities between constants k2i > kij for (i, j) 
(2, 1), fci2 > A:32, and fc23 > ki3 to obtain the simple 
answer (|58p , ([59|) , hence if we even do not use these 



^21 , ^31 ^ 
^2 ^4 



if ^31 > ^23 



, ♦ . '*^21 , ^23 , 
^23 j'^ 4 -^2 4 



«-23 ^"31 



Fig. 5. Dominant systems for case (a) (defined in Fig. |4]l 

inequalities for the further simplification, this does 
not guarantee the higher accuracy of formulas. 



5.1.3. Eigenvalues and eigenvectors 

At the next step, we have to restore and cut the 
cycles. First cycle to cut is the result of cycle gluing, 
Al ^ A3. It is necessary to delete the limiting step, 
i.e. the reaction with the smallest rate constant. If 
^31 > ^23, then we get Al ^ A3. If, inverse, ^23 > 
fcg]^, then we obtain Al ^ ^3. 

After that, we have to restore and cut the cycle 
which was glued into the vertex Al . This is the two- 
vertices cycle Al ^ A2. The limiting step for this 
cycle is Ai ^ A2, because ^21 ^ fci2. If k^i > ^23, 
then following the rule visualized by Fig. [U we get 
the dominant system Ai A2 A3 with reaction 
rate constants ^21 for Ai A2 and k^i for A2 
A3. If fc23 > k^^ then we obtain Ai — > A2 ^ A3 with 
reaction rate constants A:2i for Ai — s- A2 and /C23 for 
A2 <— A3. All the procedure is illustrated by Fig.[5l 

The eigenvalues and the correspondent eigenvec- 
tors for dominant systems in case (a) are represented 
below in zero-one asymptotic. 

(i) kl^ > k23, the dominant system Ai ^ A2 ^ 
A3, 



Ao = 0, (0,0,1), = (1,1,1) 

Ai«-fc2i, (1,-1,0), (1,0,0) 

X2^~kli, r^^ (0,1,-1), f ^{1,1,0) 

(60) 

(ii) ^23 > k^i, the dominant system Ai ^ A2 <— 
A3, 

Ao = 0, (0,1,0), = (1,1,1); 

Ai«-fc2i, (1,-1,0), (1,0,0); 

A2«-/C23, (0,-1,1), (0,0,1). 

(61) 

Here, the value of k^i is given by formula ([57|) . 
With higher accuracy, in case (a) 



w w w 
k2i ' ki2 ' A;23 



(62) 
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(b) ,.-;» ^ ^31 IT '"^13 A 

4 /ti3t^4— ►4 — ^2 

^3 if ^13 > ^31 
Fig. 6. Dominant systems for case (b) (defined in Fig. UJ 

where 




in according to ([55)) . ((55|) . 

5.2. Auxiliary system (b): A3 —> Ai ^ A2; 
ki2 > k32, ki3 > k23 

5.2.1. Gluing cycles 

The attractor is a cycle Ai ^ A2 again, and this 
is not a sink. We have to glue the cycle Ai A2 into 
one new component A\ and to add a new reaction 
A\ A3 with the rate constant k^i given by formula 
((57)) . As a result, we get an new system, Al ^ A3 
with reaction rate constants kl^ (for Al — > A3) and 
initial fci3 (for A\ ^ A3). At this stage, the only 
difference from the case (a) is the reaction Al ^ A3 
rate constant ki3 instead of fes • 

5.2.2. Steady states 

For the steady states we have to repeat formulas 
((58| (|59|) with minor changes (just use fci3 instead 

of k23). 



w = 




1 








. C3 = 




1 

^31 
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1 ■ 

fcl3 
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^21 
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1 ■ 

fel2 




ki2 



(63) 



5.2.3. Eigenvalues and eigenvectors 

The structure of the dominant system depends on 
the limiting step of the cycle Al ^ A3 (Fig. [5]). If 
k^^ > ki3 , then in the dominant system remains the 
reaction Al A3 from this cycle. After restoring 
the glued cycle Ai ^ A2 it is necessary to delete the 
slowest reaction from this cycle too. This is always 
Al <— A2, because Ai — > A2 is the fastest reaction. 
The reaction Al A3 transforms into A2 A3, 
because A2 is the head of the limiting step Ai ^ A2 
(see Fig.[T]). Finally, we get Ai A2 ^ A3. 



If fci3 > fcgj^, then then in the dominant system 
remains the reaction ^3 — > Ai, and the dominant 
system is A3 ^ Ai ^ A2 (Fig. [6]). 

The eigenvalues and the correspondent eigenvec- 
tors for dominant systems in case (b) are represented 
below in zero-one asymptotic. 

(i) kl^ > ki3, the dominant system Ai A2 ^ 
A3, 

Ao = 0, r"«(0,0,l), Z" = (l,l,l); 

Ai«-fc2i, (1,-1,0), (1,0,0); 

A2«-fc^i, r^^ (0,1,-1), (1,1,0); 

(64) 

(ii) ki3 > /cg]^, the dominant system ^.3 ^ Ai — *• 
A2, 

Ao = 0, (0,1,0), /" = (1,1,1); 

Ai«-fc2i, (1,-1,0), /I « (1,0,0); 

A2«-fci3, (0,-1,1), Z^s, (0,0,1). 

(65) 

Here, the value of fcg^ is given by formula (|57p . The 
only difference from case (a) is the rate constant ^23 
instead of ^13. 

With higher accuracy, in case (b) 

f w' w' w \ 
~ 7— , 1— , 1— , (66) 
\k2i ki2 ki3j 

where w and w' are given by formula ([55]) . 
5.3. Auxiliary system (c): Ai ^ A2 ^ A3; 

k32 > ki2, k23 > ki3 

5.3.1. Gluing cycles 

The attractor is a cycle A2 ^ A3. This is not a 
sink, because two outgoing reactions exist: A2 Ai 
and A3 ^ Al. We have to glue the cycle A2 ^ A3 
into one new component Aj and to add a new reac- 
tion A2 — > Al with the rate constant fc^j- The defini- 
tion of this new constant depends on normalized the 
steady-state distribution in this cycle. If C2, C3 are 
the steady-state concentrations (with normalization 
C2 + C3 = 1), then 

kl2 w max{fci2C2, fci3C3}. 

If we use limitation in the glued cycle explicitly, then 
we get the direct analogue of ([57)) in two versions: 
one for fc32 > A:23, another for /c23 > ^32- But we can 
skip this simplification and write 

klr^ w max{fci2W*/^32, ^13^*7^23}, (67) 



30 




A ^ if /t,, > kj 



if /t23 > 



Fig. 7. Dominant systems for case (c) (defined in Fig. |4]l 

where 

1 



1 

fc32 



1 
fe23 



5.3.2. Steady states 

Exactly as in the cases (a) and (b) we can find 
approximation of steady state using steady states in 
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(68) 



5.3.3. Eigenvalues and eigenvectors 

The hmiting step in the cycle Ai ^ A\ in known, 
this is Ai ^ A\. There are two possibilities for the 
choice on limiting step in the cycle A2 A3. If 
^32 > fes, then this limiting step is A2 ^ A3, and 
the dominant system is Ai ^ A2 ^ A-^. If fc23 > 
^32, then the dominant system is Ai ^ A2 ^ A^^ 
(Fig. El). 

The eigenvalues and the correspondent eigenvec- 
tors for dominant systems in case (b) are represented 
below in zero-one asymptotic. 

(i) ^32 > fes, the dominant system Ai A2 
A3, 



Ao = 0, r°«(0,0,l), ;° = (1,1,1) 

Ai « -fc2i, ~ (1,-1,0), « (1,0,0) 

A2«-fc32, (0,1,-1), (1,1,0) 

(69) 

(ii) ^23 > ^32, the dominant system Ai A2 ^ 

Ao = 0, r°«(0,l,0), /°-(l,l,l); 

Ai « -fc2i, ~ (1,-1,0), « (1,0,0); 

A2«-fc23, (0,-1,1), (0,0,1). 

(70) 

With higher accuracy the value of is given by 
formula of the steady state concentrations ([68]) . 



5.4. Auxiliary system (d): Ai 

k32 > ki2, fcl3 > ^23 



A. 



A3 



Ai; 



This is a simple cycle. We discussed this case in 
details several times. To get the dominant system 
it is sufficient just to delete the limiting step. Ev- 
erything is determined by the choice of the mini- 
mal constant in the couple {^32, ^13}- Formulas for 
steady state are well known too: PT|) . p^ . p^ . 

This is not necessary to discuss all orderings of 
constants, because some of them are irrelevant to 
the final answer. For example, in this case (d) in- 
terrelations between constants ^31, fc23, and ki2 are 
not important. 

5.5. Resume: zero-one multiscale asymptotic for 
the reversible reaction triangle 

We found only three topologically different ver- 
sion of dominant systems for the reversible reaction 
triangle: (i) Ai ^ A2 ^ A3, (ii) Ai ^ A2 ^ A3, 
and (iii) A3 — > Ai — > A2. Moreover, there exist only 
two versions of zero-one asymptotic for eigenvectors: 
the fastest eigenvalue is always — fc2i (because our 
choice of enumeration) , the correspondent right and 
left eigenvectors (fast mode) are: « (1,-1,0), 
l^ — (1,0,0). (The difference between systems (ii) 
and (iii) appears in the first order of the slow/fast 
constants ratio.) 

If in the steady state (almost) all mass is concen- 
trated in A2 (this means that r*' « (0, 1,0), domi- 
nant systems (ii) or (iii)), then « (0,-1,1) and 
P w (0,0,1). If in the steady state (almost) all 
mass is concentrated in A3 (this means that r*^ w 
(0, 0, 1), dominant system (i)), then w (0, 1, —1) 
and l"^ w (0,1,0). We can see that the dominant 
systems of the forms (ii) and (iii) produce the same 
zero-one asymptotic of eigenvectors. Moreover, the 
right eigenvectors « (0,1,-1) coincide for all 
cases (there is no difference between r^ and — r^), 
and the difference appears in the left eigenvector l"^. 
Of course, this peculiarity (everything is regulated 
by the steady state asymptotic) results from the sim- 
plicity of this example. 

In the zero-one asymptotic, the reversible reaction 
triangle is represented by one of the reaction mech- 
anisms, (i) or (iii). The rate constant of the first re- 
action Ai A2 is always ki2. The direction of the 
second reaction is determined by a system of linear 
uniform inequalities between logarithms of rate con- 
stants. The logarithm of effective constant of this 
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reaction is the piecewise linear function of the log- 
arithms of reaction rate constants, and the switch- 
ing between different pieces is regulated by linear 
inequalities. These inequalities are described in this 
section, and most of them are represented in Figs. ID- 
[71 One can obtain the first-order approximation of 
eigenvectors in the slow/fast constants ratio from 
the Appendix 1 formulas. 

6. Three zero-one laws and nonequilibrium 
phase transitions in multiscale systems 

6.1. Zero-one law for steady states of weakly ergodic 
reaction networks 

Let us take a weakly ergodic network W and ap- 
ply the algorithms of auxiliary systems construction 
and cycles gluing. As a result we obtain an auxiliary 
dynamic system with one fixed point (there may be 
only one minimal sink). In the algorithm of steady 
state reconstruction (Subsec. l43)) we always operate 
with one cycle (and with small auxiliary cycles in- 
side that one, as in a simple example in Subsec. [2T9|) . 
In a cycle with limitation almost all concentration 
is accumulated at the start of the limiting step , 
(fT4)) . Hence, in the whole network almost all concen- 
tration will be accumulated in one component. The 
dominant system for a weekly ergodic network is an 
acyclic network with minimal element. The minimal 
element is such a component A„iin that there exists 
a oriented path in the dominant system from any 
element to ^min- Almost all concentration in the 
steady state of the network W will be concentrated 
in the component Amin- 

6.2. Zero-one law for nonergodic multiscale 
networks 

The simplest example of nonergodic but con- 
nected reaction network is <— ^2 ^ ^3 with 
reaction rate constants fci , /c2 ■ For this network, in 
addition to 5°(c) — ci -\- C2 -\- 03 a. kinetic conserva- 
tion law exist, b^(c) ~ §^ — p.. xhe result of time 
fci 

evolution, lim^^^oo exp(ii't)c (|30p . is described by 
simple formula ([?T|) : 

lim exp(Xt)c = b\c){l, 0, 0) -I- &^(c)(0, 1, 1), 

t — *oo 

where b^c) + b^{c) = b"{c) and ^^^b^c) - 
^^b^{c) = b^ic). If fci > k2 then b^{c) ^ ci + C2 
and &^(c) w C3. If fci ^ fc2 then b^{c) ~ Ci and 



6^(c) « C2 + C3. This simple zero-one law (either 
almost all amount of A2 transforms into Ai , or al- 
most all amount of A2 transforms into A3) can be 
generalized onto all nonergodic multiscale systems. 

Let us take a multiscale network and perform the 
iterative process of auxiliary dynamic systems con- 
struction and cycle gluing, as it is prescribed in Sub- 
sec. 14.31 After the final step the algorithm gives the 
discrete dynamical system (p™ with fixed points AJ^. 

The fixed points AJ\ of the discrete dynamical sys- 
tem are the glued ergodic components Gi C A 
of the initial network W. At the same time, these 
points are attractors of $™. Let us consider the cor- 
respondent decomposition of this system with parti- 
tion A™' — UiAtt{A'jl). In the cycle restoration dur- 
ing construction of dominant system dommod(yV) 
this partition transforms into partition of ^ = 
UiUi, Att{AJ\) transforms into Ui and Gi C Ui (and 
Ui transforms into Att{AJ\) in hierarchical gluing of 
cycles). 

It is straightforward to see that during construc- 
tion of dominant systems for W from the network 
V™ no connection between Ui are created. There- 
fore, the reaction network dom mod(W) is a union of 
networks on sets Ui without any link between sets. 

If Gi , . . . Gm are all ergodic components of the 
system, then there exist m independent positive lin- 
ear functionals 6^(c), ... fe™(c) that describe asymp- 
totical behaviour of kinetic system when t —^ 00 
([30)) . For dommod(yV) these functionals are: fe'(c) = 
'^AeUi ^A where ca is concentration of A. Hence, for 
the initial reaction network W with well separated 
constants 

b\c) « (71) 

Aeu, 

This is the zero-one law for multiscale networks: for 
any /, i, the value of functional 6' ([5(1)) on basis vector 
e\ &'(e*), is either close to one or close to zero (with 
probability close to 1). We already mentioned this 
law in discussion of a simple example ()3ip . The ap- 
proximate equality ()7ip means that for each reagent 
A G A there exists such an ergodic component G of 
W that A transforms when i — > 00 preferably into 
elements of G even if there exist paths from A to 
other ergodic components of W. 

6.3. Dynamic limitation and ergodicity boundary 

Dominant systems are acyclic. All the stationary 
rates in the first order are limited by limiting steps 
of some cycles. Those cycles are glued in the hier- 
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archical cycle gluing procedure, and their limiting 
steps are deleted in the cycles surgery procedures 
(see Subsec.lOland Fig. [J). 

Relaxation to steady state of the network is multi- 
exponential, and now we are interested in estimate 
of the longest relaxation time r: 



r = l/mm{-ReXi\X^ ^ 0} 



(72) 



Is there a constant that limits the relaxation time? 
The general answer for multiscale system is; 1/r is 
equal to the minimal reaction rate constant of the 
dominant system. It is impossible to guess a priori, 
before construction of the dominant system, which 
constant it is. Moreover, this may be not a rate con- 
stant for a reaction from the initial network, but a 
monomial of such constants. 

Nevertheless, sometimes it is possible to point the 
reaction rate constant that is limiting for the re- 
laxation in the following sense. For known topology 
of reaction network and given ordering of reaction 
rate constants we find such a constant (ergodicity 
boundary) that 

with a < 1 is a function of constants kj > k^-- This 
means that l/k-r gives the lower estimate of the re- 
laxation time, but r could be larger. In addition, we 
show that there is a zero-one alternative too: if the 
constants are well separated then either a ~ 1 or 
a < 1. 

We study a multiscale system with a given reac- 
tion rate constants ordering, fej^ > kj^ > . . . > kj^^ . 
Let us suppose that the network is weakly ergodic 
(when there are several ergodic components, each 
one has its longest relaxation time that can be found 
independently) . We say that kj^ , 1 < r < n is the er- 
godicity boundary kr if the network of reactions with 
parameters kj-^ , /cj^ , . . . , kj^ (when fcj^^^ = ■■■kj^ = 
0) is weakly ergodic, but the network with parame- 
ters fcjj ,kj2,..., kj^-i (when kj^ — kj^_^-^ = ...kj^ = 
0) it is not. In other words, when eliminating re- 
actions in decreasing order of their characteristic 
times, starting with the slowest one, the ergodicity 
boundary is the constant of the first reaction whose 
elimination breaks the ergodicity of the reaction di- 
graph. This reaction we also call the "ergodicity 
boundary" . 

Let us describe the possible location of the er- 
godicity boundary in the general multiscale reac- 
tion network (W) . After deletion of reactions with 
constants kj^ , kj^^-^ , ...fcj„ from the network two er- 



. ^ ^3 A2 A3 

a) b) 



Fig. 8. Two basic examples of ergodicity boundary reaction: 
(a) Connection between ergodic components; (b) Connection 
from one ergodic component to element that is connected 
to the both ergodic components by oriented paths. In both 
cases, for £ = 0, the ergodic components are {A2} and {A3}. 

godic components (minimal sinks) appear, Gi and 
G2 ■ The ergodicity boundary starts in one of the er- 
godic components, say Gi, and ends at the such a 
reagent B that another ergodic component, G2, is 
reachable by B (there exists an oriented path from 
B to some element of G2). 

An estimate of the longest relaxation time can 
be obtained by applying the perturbation theory 
for linear operators to the degenerated case of the 
zero eigenvalue of the matrix K. We have K = 



K^r {kji , kj 



, kj^-i) -\- kj^Q + o{kr), where K^r 



is obtained from K by letting kr = fc^+i — . . .kn = 
0, Q is a constant matrix of rank 1, and o{kr) in- 
cludes terms that are negligible relative to kr- The 
zero eigenvalue is twice degenerate in K^r and not 
degenerate in K^r + krQ. One gets the following es- 
timate: 

a^>r>a^, (74) 

where a, a > are some positive functions of 
ki, k2, ■ ■ ■ , kr-i (and of the reaction graph topol- 
ogy). 

Two simplest examples demonstrate two types of 
dependencies of t on kr'. 

(i) For the reaction mechanism Fig.[S^ 

min{— i?eA| — e, 

\i e < ki + k2- 

(ii) For the reaction mechanism Fig.[5j3 

min{— i?eA} = efc2/(fci + ^2) + o{e), 

if £ < fci + /e2. For well separated parameters 
there exists a zero-one (trigger) alternative: if 
ki <g; ^2 then miTi\^Q{—ReX\ « S] if, inverse, 
ki ^ ^2 then TmTi\^Q{— ReX\ = o(e). 
In general multiscale network, two type of obsta- 
cles can violate approximate equality r w l/fcr- Fol- 
lowing the zero-one law for nonergodic multiscale 
networks (previous subsection) we can split the set 
of all vertices into two subsets, Ui and U2- The dom- 
inant reaction network dommod(yV) is a union of 
networks on sets Ui 2 without any link between sets. 



33 



If the ergodicity boundary reaction starts in the 
ergodic component Gi and ends at B which does 
not belong to the "opposite" basin of attraction U2, 
then T 3> 1/kr- This is the first possible obstacle. 

Let the ergodicity boundary reaction start at A e 
Gi and end at i? G U2- We define the maximal linear 
chain of reactions in dominant system with start at 
B: B ^ ... . This chain belongs to C/2- Let us extend 
this chain from the left by the ergodicity bound- 
ary: A ^ B ^ ... . Relaxation time for the net- 
work of r reactions (with the kinetic matrix K<:r = 



K^rikj^ , kj. 



kj^Q) is, approximately, 



the relaxation time of this chain, i.e. where k 
is the minimal constant in the chain. There may ap- 
pear a monomial constant fc <C fci-. In that case, 
T 3> 1/^7- , and relaxation is limited by this mini- 
mal k or by some of constants kj^ , p > r or by some 
of their combinations. This existence of a monomial 
constant fc ^ fc,- in the maximal chain A B 
... from the dominant system is the second possible 
obstacle for approximate equality r « l/k-r-. 

If there is neither the first obstacle, nor the second 
one, then r « 1/kr. The possibility of these obsta- 
cles depends on the definition of multiscale ensem- 
bles we use. For example for the log-uniform distri- 
bution of rate constants in the ordering cone kj-^ > 
kj^ > ... > kj^ (Sec. |3|) the both obstacles have 
nonzero probability, if they are topologically possi- 
ble. On the other hand, if we study asymptotic of re- 
laxation time at e — > for ki^, — ekj^^i for given val- 



ues of kj^ , kj2 , . 



1 k^. 



then for sufficiently small 



e > the second obstacle is impossible. 

Thus, the well known concept of stationary reac- 
tion rates limitation by "narrow places" or "limiting 
steps" (slowest reaction) should be complemented 
by the ergodicity boundary limitation of relaxation 
time. It should be stressed that the relaxation pro- 
cess is limited not by the classical limiting steps (nar- 
row places), but by reactions that may be absolutely 
different. The simplest example of this kind is an ir- 
reversible catalytic cycle: the stationary rate is lim- 
ited by the slowest reaction (the smallest constant), 
but the relaxation time is limited by the reaction 
constant with the second lowest value (in order to 
break the weak ergodicity of a cycle two reactions 
must be eliminated). 



6.4. Zero-one law for relaxation modes 
( eigenvectors ) and lumping analysis 

For kinetic systems with well-separated constants 
the left and right eigenvectors can be explicitly es- 
timated. Their coordinates are close to ±1 or 0. We 
analyzed these estimates first for linear chains and 
cycles ([IJ and then for general acyclic auxiliary dy- 
namical systems ^ ([351), (EZl- The distri- 
bution of zeros and ±1 in the eigenvectors compo- 
nents depends on the rate constant ordering and 
may be rather surprising. Perhaps, the simplest ex- 
ample gives the asymptotic equivalence (for k~ ^ 
/ci, A:i_|_i) of the reaction network <-> Ai+i — s- Ai+2 
with rate constants ki, k~ and fc^+i to the reaction 
network Aij^i ^ — s- with rate constants k~ 
(for the reaction Ai+i — > Ai) and ki+iki/k~ (for the 
reaction Ai — > ^i+2) presented in in Subsec. [2Jl 

For reaction networks with well-separated con- 
stants coordinates of left eigenvectors T are close 
to or 1. We can use the left eigenvectors for co- 
ordinate change. For the new coordinates Zi — Pc 
(eigenmodes) the simplest equations hold: Zi = XiZi. 
The zero-one law for left eigenvectors means that the 
eigenmodes are (almost) sums of some components: 
Zi = X^jgy '^i f*^^ some sets of numbers Vi. Many ex- 
amples, dll, ([38|) . (|5T|l . ([54|l . demonstrate that some 
of Zi can include the same concentrations: it may be 
that Vi nVj ^ for some i ^ j. Aggregation of 
some components (possibly with some coefficients) 
into new group components for simplification of ki- 
netics is the major task of lumping analysis. 

Wei and Kuo studied conditions for exact 
(IWei fc Kuol (|l969f )) and approximate (|Kuo fc Weil 
(Il969l )) linear lumping. More recently, sensitiv- 
ity analysis and Lie group approach w e re a 
plied to lumping analysis I Li fc Rabita (|l98' 
Toth. Li. Rabitz. k TomlinI (|l997^ )■ and more 
general nonlinear forms of lumped concentra- 
tions are used (for example, Zi could be ratio- 
nal function of c). The power of lumping us- 
ing a tim e-scale based approach was demon - 
strated by IWhitehouse Tomlin. fc Pillinel (|2004l ): 
Liao fc Lightfooll ( 19881 ). This computationally 
cheap approach combines ideas of sensitivity analy- 
sis with simple and useful grouping of species with 
similar lifetimes and similar topological properties 
caused by connections of the species in the reaction 
networks. The lumped concentrations in this ap- 
proach are simply sums of concentrations in groups. 

Kinetics of multiscale systems studied in this 
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paper and developed theory of dynamic limitation 
demonstrates that in multiscale limit lumping anal- 
ysis can work (almost) exactly. Lumped concentra- 
tions are sums in groups, but these groups can in- 
tersect and usually there exist several intersections. 

6.5. Nonequilibrium phase transitions in multiscale 
systems 

For each zero-one law specific sharp transitions 
exist: if two systems in a one-parametric family have 
different zero-one steady states or relaxation modes, 
then somewhere between a point of jump exists. Of 
course, for given finite values of parameters this will 
be not a point of discontinuity, but rather a thin zone 
of fast change. At such a point the dominant sys- 
tem changes. We can call this change a nonequilib- 
rium phase transition. Here we identify a "multiscale 
nonequilibrium phase" with a dominant system. 

A point of phase transition can be a point where 
the order of parameters changes. But not every 
change of order causes the change of dominant 
systems. On the other hand, change of order of 
some monomials can change the dominant system 
even if the order of parameters persists (examples 
are presented in previous section). Evolution of a 
parameter-dependent multiscale reaction network 
can be represented as a sequence of sharp change 
of dominant system. Between such sharp changes 
there are periods of evolution of dominant system 
parameters without qualitative changes. 

7. Limitation in modular structure and 
solvable modules 

7.1. Modular limitation 



modular structure of reaction network. Let for the 
reaction network W the set of elementary reactions 
TZ is partitioned on some modules: TZ = UiTZi. We 
can consider the related multiscale ensemble of re- 
action constants: let the ratio of any two rate con- 
stants inside each module be bounded (and sepa- 
rated from zero, of course), but the ratios between 
modules form a well separated ensemble. This can 
be formalized by multiplication of rate constants of 
each module TZi on a time scale coefficient fc^. If we 
assume that In ki are uniformly and independently 
distributed on a real line (or ki are independently 
and log- uniformly distributed on a sufficiently large 
interval) then we come to the problem of modular 
limitation. The problem is quite general: describe 
the typical behavior of multiscale ensembles for sys- 
tems with given modular structure: each module has 
its own time scale and these time scales are well sep- 
arated. 

Development of such a general theory is outside 
the scope of our paper, and here we just find building 
blocks for the future theory, solvable reaction mod- 
ules. There may be many various criteria of selec- 
tion the reaction modules. Here are several possible 
choices: individual reactions (we developed the the- 
ory of multiscale ensembles of individual reactions 
in this paper), couples of mutually inverse reactions, 
as we mentioned above, acyclic reaction networks. 

Among the possible reasons for selection the class 
of reaction mechanisms for this purpose, there is one 
formal, but important: the possibility to solve the 
kinetic equation for every module in explicit analyt- 
ical (algebraic) form with quadratures. We call these 
systems "solvable". 

7.2. Solvable reaction mechanisms 



The simplest one-constant limitation concept can- 
not be applied to all systems. There is another very 
simple case based on exclusion of "fast equilibria" 
Ai ^ Aj. In this limit, the ratio of reaction con- 
stants Kij — kij /kji is bounded, < a < Kij < b < 
oo, but for different pairs {i,j), {I, s) one of the in- 
equalities kij <^ kis or kij 3> kis holds. (One usually 
calls these K "equilibrium con stant" , even if there is 
no relevant thermodynamics, discussed 
that case systematically for some real examples. Of 
course, it is possible to create the theory for that 
case very similarly to the theory presented above. 
This should be done, but it is worth to mention now 
that the limitation concept can be applied to any 



Let us describe all solvable reaction systems (with 
mass action law), linear and nonlinear. 

Formally, we call the set of reaction solvable, if 
there exists a linear transformation of coordinates 
c !■ a such that kinetic equation in new coordinates 
for all values of reaction constants has the triangle 
form: ^ 

-rr = ft{ai,a2,...ai). (75) 
at 

This system has the lower triangle Jacobian matrix 
dcLi I doj . 

To construct the general mass action law system 
we need: the list of components, A = {Ai, ...An} 
and the list of reactions (the reaction mechanism) : 
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E 



(76) 



where r is the reaction number, a^i, /3rfe are nonneg- 
ative integers (stoichiometric coefficients). Formally, 
it is possible that all f3k — or all Oi ~ 0. We allow 
such reactions. They can appear in reduced models 
or in auxiliary systems. 

A real variable Ci is assigned to every component 
Ai^ Ci is the concentration ol Ai, c is the concentra- 
tion vector with coordinates Cj. The reaction kinetic 
equations are 



(77) 



where "fr is the reaction stoichiometric vector with 
coordinates jri = Pri — C(ri, Wr{c) is the reaction 
rate. For mass action law, 



n (c) — kr JJ^ ( 



(78) 



where kr is the reaction constant. 

Physically, equations ([77)) correspond to reactions 
in fixed volume, and in more general case a multiplier 
V (volume) is necessary: 



A{Vc) 
dt 



Here we study the systems ([77]) and postpone any 
further generalization. 

The first example of solvable systems give the sets 
of reactions of the form 



cXriAi 



(79) 



(components Ak on the right hand side have higher 
numbers k than the component Ai on the left hand 
side, i < k). For these systems, kinetic equations 
()77p have the triangle form from the very beginning. 

The second standard example gives the couple of 
mutually inverse reactions: 



fc, 



(80) 



these reactions have stoichiometric vectors ±7, 7^ = 
Pi — ai. The kinetic equation c = (w^ — w^)^ has the 
triangle form (|75p in any orthogonal coordinate sys- 
tem with the last coordinate an = (7, c) = 7iCj. 
Of course, if there are several reactions with propor- 
tional stoichiometric vectors, the kinetic equations 
have the triangle form in the same coordinate sys- 
tems. 



The general case of solvable systems is essen- 
tially a combination of that two (|79|) . ([80| . with 
some generalization. Here we follow the book by 
Gorban. Bvkov. fc Yablonskiil ( 19861 ) and present 



an algorithm for analysis of reaction network solv- 
ability. First, we introduce a relation between re- 
actions "rth reaction directly affects the rate of 
sth reaction" with notation r ^ s: r — > s if there 
exists such Ai that ^riOisi 7^ 0. This means that 
concentration of Ai changes in the rth reaction 
{iri 7^ 0) and the rate of the sth reaction depends 
on Ai concentration [asi 7^ 0) . For that relation we 
use r — !■ s. For transitive closure of this relation we 
use notation r > s { "rth reaction affects the rate of 
sth reaction" ) : r ^ s if there exists such a sequence 

Sl, 32, ... Sq that r Si — > S2 — > ... Sq s. 

The hanging component of the reaction network 
yy is such Ai £ A that for all reactions a^i = 0. This 
means that all reaction rates do not depend on con- 
centration of Ai. The hanging reaction is such ele- 
ment of TZ with number r that r )^ s only if 7s = X'jr 
for some number A. An example of hanging compo- 
nents gives the last component An for the triangle 
network (j79|) . An example of hanging reactions gives 
a couple of reactions ([80]) if they do not affect any 
other reaction. 

In order to check solvability of the reaction net- 
work W we should find all hanging components and 
reactions and delete them from A and TZ, corre- 
spondingly. After that, we get a new system, Wi 
with the component set Ai and the reaction set TZi. 
Next, we should find all hanging components and 
reactions for Wi and delete them from Ai and TZi. 
Iterate until no hanging components or hanging re- 
actions could be found. If the final set of components 
is empty, then the reaction network W is solvable. 
If it is not empty, then W is not solvable. 

For example, let us consider the reaction mecha- 
nism with A — {Ai, ^2, ^3, A4} and reactions Ai + 
A2 2A3, A1+A2 As + Ai, A3 Ai, Ai A3. 
There are no hanging components, but two hanging 
reactions, A3 A4 and A4 A3. After deletion 
of these two reactions, two hanging components ap- 
pear, ^3 and A4. After deletion these two compo- 
nents, we get two hanging reactions, Ai + A2 0, 
Ai + A2 ^ (they coincide) . We delete these reac- 
tions and get two components Ai, A2 without reac- 
tions. After deletion these hanging components we 
obtain the empty system. The reaction network is 
solvable. 

An oriented cycle of the length more than two is 
not solvable. For each number of vertices one can cal- 
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(i) M 

(ii) A, 



culate the set of all maximal solvable mechanisms. 
For example, for five components there are two max- 
imal solvable mechanisms of monomolecular reac- 
tions: 

A2 Ai, Ai ^4, A2 ^3, Ai -> 
A, Ai ^ A5, Ai ^ A5; 
A2, Ai A3, Ai -> Ai, Ai A5, 
A2 A3, Ai ^ A5. 
It is straightforward to check solvability of these 
mechanism. The first mechanism has a couple of 
hanging reactions, Ai ^ A^. After deletion of 
these reactions, the system becomes acyclic, of the 
form (j79p . The second mechanism has two couples 
of hanging reactions, A2 A3 and Ai ^ A^. 
After deletion of these reactions, the system also 
transforms into the triangle form ((79|) . It is im- 
possible to add any new monomolecular reactions 
between {Ai, A2, A3, Ai, A^} to these mechanisms 
with preservation of solvability, and any solvable 
monomolecular reaction network with five reagents 
is a subset of one of these mechanisms. 

Finally, we should mention connections between 
solvable reaction n etworks and so lvable Lie algebras 
(|jacobsonl (|l979t ): lde Graaj (f2000,l ). Let us remind 



that matrices Mi, ... Mq generate a solvable Lie al- 
gebra if and only if they could be transformed simul- 
taneously into a triangle form by a change of basis. 

The Jacobian matrix for the mass action law ki- 
netic equation ((77l) is: 



^)=T.-rJr = Y.^M.,, (81) 



where 



Jr = 7rajdiag | 



Mrj = Urj^rG^^ , 



1 1 1 

Cl ' C2 ' C„ 



(82) 



aj is the vector row {ari, ■■■ am), e-'^ is the jth 
basis vector row with coordinates e|, = Sjk ■ 

The Jacobian matrix (f8T|) should have the lower 
triangle form in coordinates (|75p for all nonneg- 
ative values of rate constants and concentrations. 
This is equivalent to the lower triangle form of all 
matrices Mrj in these coordinates. Because usually 
there are many zero matrices among Mrj , it is con- 
venient to describe the set of nonzero matrices. 

For the rth reaction = {i \ art ^ 0}. The reac- 
tion rate Wr depends on ci if and only if z G J^. For 
each z = 1 , . . . n we define a matrix 



0,0,. 



7r 



The ith column of this matrix coincides with the vec- 
tor column 7r. Other columns are equal to zero. For 
each r we define a set of matrices M.r = {mri \ i € 
Ir}, and Ai = UrMr- The reaction network W is 
solvable if and only if the finite set of matrices M. 
generates a solvable Lie algebra. 

Classification of finite dimension al solvable Lie al- 
gebra s remains a difficult problem (|de Graai (|2000l 
20051 )'!. It seems plausible that the classification of 
solvable algebras associated with reaction networks 
can bring new ideas into this field of algebra. 

8. Conclusion: Concept of limit 
simplification in multiscale systems 

In this paper, we study networks of linear reac- 
tions. For any ordering of reaction rate constants we 
look for the dominant kinetic system. The dominant 
system is, by definition, the system that gives us the 
main asymptotic terms of the stationary state and 
relaxation in the limit for well separated rate con- 
stants. In this limit any two constants are connected 
by the relation ^ or 

The topology of dominant systems is rather sim- 
ple; they are those networks which are graphs of dis- 
crete dynamical systems on the set of vertices. In 
such graphs each vertex has no more than one out- 
going reaction. This allows us to construct the ex- 
plicit asymptotics of eigenvectors and eigenvalues. 
In the limit of well separated constants, the coordi- 
nates of eigenvectors for dominant systems can take 
only three values: ±1 or 0. All algorithms are repre- 
sented topologically by transformation of the graph 
of reaction (labeled by reaction rate constants) . We 
call these transformations "cycles surgery" , because 
the main operations are gluing cycles and cutting 
cycles in graphs of auxiliary discrete dynamical sys- 
tems. 

In the simplest case, the dominant system is de- 
termined by the ordering of constants. But for suffi- 
ciently complex systems we need to introduce aux- 
iliary elementary reactions. They appear after cycle 
gluing and have monomial rate constants of the form 
kt; = Yii ^i* ■ The dominant system depends on the 
place of these monomial values among the ordered 
constants. 

Construction of the dominant system clarifies 
the notion of limiting steps for relaxation. There is 
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an exponential relaxation process that lasts much 
longer than the others in (|44| , ((53| . This is the slow- 
est relaxation and it is controlled by one reaction in 
the dominant system, the limiting step. The limit- 
ing step for relaxation is not the slowest reaction, or 
the second slowest reaction of the whole network, 
but the slowest reaction of the dominant system. 
That limiting step constant is not necessarily a re- 
action rate constant for the initial system, but can 
be represented by a monomial of such constants as 
well. 

The idea of dominant subsystems in asymptotic 
analysis was pr oposed by Newton and developed by 
Kruskall (jl963r ). A modern introduction with some 
historical review is presented by IWhitd (,20061 In 
our analysis we do not use the powers of small pa- 
rameters (as it was d o ne bv IVishik fc Liusternik 



1960): iLidskr 



20041 ): iKruska 




.Akian. Bapat. fc Gaubert 
Whitel (120061 )1. but operate 



directly with the rate constants ordering. 

To develop the idea of systems with well separated 
constants to the state of a mathematical notion, we 
introduce multiscale ensembles of constant tuples. 
This notion allows us to discuss rigorously uniform 
distributions on infinite space and gives the answers 
to a question: what does it mean "to pick a multi- 
scale system at random" . 

Some of results obtained are rather surprising and 
unexpected. First of all is the zero-one asymptotic 
of eigenvectors. Then, the good approximation to 
eigenvectors does not give approximate eigenvectors 
(the inverse situation is more common: an approx- 
imate eigenvector could be far from the eigenvec- 
tor). The almost exact lumping analysis provided 
by the zero-one approximation of eigenvectors has 
an unexpected property: the lumped groups for dif- 
ferent eigenvalues can intersect. Rather unexpected 
seems the change of reaction sequence when we con- 
struct the dominant systems. For example, asymp- 
totic equivalence (for 3> fcj, /c^+i) of the reaction 
network Ai <-> ^i+i — > ^^+2 with rate constants 
ki, k~ and fc^+i to the reaction network Ai^i — > 
Ai — > Ai^2 with rate constants k^ (for the reaction 
Ai^i — > Ai) and ki^iki/k'' (for the reaction Ai — s- 
Ai+2) is simple, but surprising (Subsec. l2.9p . And, 
of course, it was surprising to observe how the dy- 
namics of linear multiscale networks transforms into 
the dynamics on finite sets of reagent names. 

Now we have the complete theory and the exhaus- 
tive construction of algorithms for linear reaction 
networks with well separated rate constants. There 
are several ways of using the developed theory and 



algorithms: 

(i) For direct computation of steady states and re- 
laxation dynamics; this may be useful for com- 
plex systems because of the simplicity of the 
algorithm and resulting formulas and because 
often we do not know the rate constants for 
complex networks, and kinetics that is ruled 
by orderings rather than by exact values of 
rate constants may be very useful; 

(ii) For planning of experiments and mining the 
experimental data - the observable kinetics is 
more sensitive to reactions from the dominant 
network, and much less sensitive to other re- 
actions, the relaxation spectrum of the domi- 
nant network is explicitly connected with the 
correspondent reaction rate constants, and 
the eigenvectors ("modes") are sensitive to 
the constant ordering, but not to exact values; 

(iii) The steady states and dynamics of the dom- 
inant system could serve as a robust first ap- 
proximation in perturbation theory or as a 
preconditioning in numerical methods. 

The developed methods are computationally 
cheap, for example, the algorithm for construction 
of dominant system has linear complexity num- 
ber of reactions). From a practical point of view, 
it is attractive to use exact rational expressions 
for the dominant system modes ([3]), ([M]) . ([55)1 in- 
stead of the zero-one approximation. Also, we can 
use exact formula (fTTj) for irreversible cycle steady 
state instead of linear approximation (fT3|) . These 
improvements are computationally cheap and may 
enhance accuracy of computations. 

From a theoretical point of view the outlook is 
more important. Let us answer the question: what 
has to be done, but is not done yet? Three directions 
for further development are clear now: 

(i) Construction of dominant systems for the re- 
action network that has a group of constants 
with comparable values (without relations ^ 
between them) . We considered cycles with sev- 
eral comparable constants in Sec. [21 but the 
general theory still has to be developed. 

(ii) Construction of dominant systems for reac- 
tion networks with modular structure. We can 
assume that the ratio of any two rate con- 
stants inside each module be bounded and sep- 
arated from zero, but the ratios between mod- 
ules form a well separated ensemble. A reac- 
tion network that has a group of constants 
with comparable values gives us an example 
of the simplest modular structure: one module 
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includes several reactions and other modules 
arise from one reaction. In Sec. [7] we describe 
all solvable modules such that it is possible 
to solve the kinetic equation for every mod- 
ule in explicit analytical (algebraic) form with 
quadratures (even for nonconstant in time re- 
action rate constants), 
(iii) Construction of dominant systems for nonlin- 
ear reaction networks. The first idea here is 
the representation of a nonlinear reaction as a 
pseudomonomolecular reaction: if for reaction 
A + B ... concentrations ca and cb are well 
separated, say, ca 2> cb, then we can consider 
this reaction as -B — s- ... with rate constant 
dependent on ca- The relative change of ca 
is slow, and we can consider this reaction 
as pseudomonomolecular until the relation 
ca ^ Cb changes to ca ^ cb- We can assume 
that in the general case only for small fraction 
of nonlinear reactions the pseudomonomolec- 
ular approach is not applicable, and this set 
of genuinely nonlinear reactions changes in 
time, but remains small. For nonlinear sys- 
tems, even the realization of the limiting step 
idea for steady states of a one-route mecha- 
nism of a catalytic reaction is nontrivial and 
was develop ed through the concep t of kin etic 
polynomial ( Lazman fc Yablonskii ( 19881 ) V 
Finally, the concept of "limit simplification" will 
be developed. For multiscale nonlinear reaction net- 
works the expected dynamical behaviour is to be 
approximated by the system of dominant networks. 
These networks may change in time but remain 
small enough. 

This hypothetical picture should give an answer 
to a very practical question: how to describe ki- 
netics beyond the standard quasi- steady-state and 
quasiequilibrium approximations ( Schnel l fc Mainj 
(|2002) ). We guess that the answer has the follow- 
ing form: during almost all time almost everything 
could be simplified and the whole system behaves 
as a small one. But this picture is also nonstation- 
ary: this small system change in time. Almost always 
"something is very small and something is very big" , 
but due to nonlinearity this ordering can change 
in time. The whole system walks along small sub- 
systems, and constants of these small subsystems 
change in time under control of the whole system 
state. The dynamics of this walk supplements the 
dynamics of individual small subsystems. 

The corresponding structure of fast-slow time 
separation in phase space is not necessarily a smooth 



slow invariant manifold, but may be similar to a 
"crazy quilt" and may consist of fragments of var- 
ious dimensions that do not join smoothly or even 
continuously. 



Appendix 1: Estimates of eigenvectors for diago- 
nally dominant matrices with diagonal gap condi- 
tion 

The famous Gershgorin theorem gives estimates 
of eigenvalues. The estimates of correspondent 
eigenvectors are not so well known. In the paper 
we use some estimates of eigenvectors of kinetic 
matrices. Here we formulate and prove these esti- 
mates for general matrices. Below A = (a^) is a 
complex n X n matrix. Pi = X^jj^^i l^-u l (sums of 
non-diagonal elements in rows), Qi = J2j j^i 
(sums of non-diagonal ele ments in columns ). 

Gershgorin theorem ( Marcus fc Mind (|l992 ). 
p. 146): The characteristic roots of A lie in the 
closed region of the z-plane 



= U (Gf = {z\\z- au\ < m. (83) 

i 

Analogously, the characteristic roots of A lie in the 
closed region G"^ of the z-plane 



G^ = [jGf iGf = { 



z \ \z — a. 



^\<Q^}■ (84) 



Areas Gf and Gf are the Gershgorin discs. 

Gershgorin discs Gf {i = 1, . . .n) are isolated, if 
Gf n Gf = for i ^ j. If discs Gf {i ^ l,...n) 
are isolated, then the spectrum of A is simple, and 
each Gershgorin disc Gf contains o ne an d only one 
eigenvalue of A (jMarcus fc Mind (|l992l ). p. 147). 
The same is true for discs Gf . 

Below we assume that Gershgorin discs Gf [i = 
1, . . . n) are isolated, this means that for all i, j 



Let us introduce the following notations: 



^33 I 



Xij 



(85) 



(86) 



mm j 

3 \aii\ 



Usually, we consider Si and Xij as sufficiently small 
numbers. In contrary, gi should not be small, (this 
is the gap condition). For example, if for any two 
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diagonal elements a^i, ajj either an 3> ajj or an <^ 
ajj, then ffi ^ 1 for all i. 

Let Ai € be the eigenvalue of A (|Ai — an] < 
Qi). Let us estimate the correspondent right eigen- 
vector x'-^' = {xi): Ax'^^'' — Aia;(^^ We take xi — \ 
and write equations for Xi (17^!): 

(an ~ ail - 9i)xi + ^ aijXj = -an, (87) 

where = Ai — an, |0i| < Qi. 
Let us introduce new variables 

x^{xi), Xi ^ Xi{aii - an) {i = 2,...n). 

In these variables, 



1 - 



an — flu 



Ea. 
n - 



ajj — On 



-Oil, 



or in matrix notations: (1 — B)x = — ai, where hi 
is a vector column with coordinates an. because of 
gap condition and smallness of Si and Xij Xij 
can consider matrix _B as a small matrix, for assume 
that ||B|| < 1 and (1 — B) is reversible (for detailed 
estimate of see below). 
For x we obtain: 



i = -hi - B{\- B)-^hi, 
and for residual estimate 



\B{1 - B)-'hi\\ < 



\B\\ 



ail 



1-P„ 

For eigenvector coordinates we get from 
an {B{\-B)-^hi\ 



an — an an ~ an 

and for residual estimate 

mi-B)-^hi).\ ^ \\B\\ \\hi\\ 



\aii - an 



1 - \au - an 



(89) 
(90) 

(91) 
(92) 



Let us give more detailed estimate of residual. For 
vectors we use h norm: ||x|| ~ The corre- 

spondent operator norm of matrix B is 

||_B|| = max < > max|6ij|. 



With the last estimate for matrix B ([55]) we find: 

Qi 



\bn\ < I I <-<-, 

\aii-an\ gi g 



(93) 



\a]]-aii\ 9j g 

where e = max^e^, x = maxij Xij-, 9 ^ min^gi. By 
definition, e > Xi a-nd for all i, j the simple estimate 



holds: \bij\ < e/g. Therefore, ||-Ba;|| < ne/g and, 
||B||/(1 — ||i?||) < ne/{g — ne) (under condition g > 
ne). Finally, ||ai|| = Qi and for residual estimate we 
get: 

-2 



an 



< 



ne 



(94) 



an gig - ne) 

More accurate estimate can be produced from in- 
equalities , if it is necessary. For our goals it is 
sufficient to use the following consequence of 

.2 



ne 



{^^l). 



(95) 



g g{g - ne) 

With this accuracy, eigenvectors of A coincide with 
standard basis vectors, i.e. with eigenvectors of di- 
agonal part of A, diagjan, . . . an„}. 



Appendix 2: Time separation and averaging in 
cycles 

In Sec. [2] we analyzed relaxation of a simple cycle 
with limitation as a perturbation of the linear chain 
relaxation by one more step that closes the chain 
into the cycle. The reaction rate constant for this 
perturbation is the smallest one. For this analysis we 
used explicit estimates ([5]) of the chain eginvectors 
for reactions with well separated constants. 

Of course, one can use estimates (p4)) . ((35|) (p6| 
and (|37[) to obtain a similar perturbation analysis 
for more general acyclic systems (instead of a linear 
chain). If we add a reaction to an acyclic system 
(after that a cycle may appear) and assume that 
the reaction rate constant for additional reaction is 
smaller than all other reaction constants, then the 
generalization is easy. 

This smallness with respect to all constants is re- 
quired only in a very special case when the addi- 
tional reaction has a form Ai — > Aj (with the rate 
constant kji) and there is no reaction of the form 
Ai ... in the non-perturbed system. In Sec. [7] 
and Appendix 1 we demonstrated that if in a non- 
perturbed acyclic system there exists another reac- 
tion of the form Ai — > ... with rate constant Ki, then 
we need inequality kji <ti Ki only. This inequality al- 
lows us to get the uniform estimates of eigenvectors 
for all possible values of other rate constants (under 
the diagonally gap condition in the non-perturbed 
system). 

For substantiation of cycles surgery we need addi- 
tional perturbation analysis for zero eigenvalues. Let 
us consider a simple cycle Ai ^ A2 ^ ... ^ An 
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Ai with reaction Ai — > ... rate constants k^. We add 
a perturbation Ai (from Ai to nothing) with 
rate constant eKi. Our goal is to demonstrate that 
the zero eigenvalue moves under this perturbation 
to Ao = —ew*{l + Xw), the correspondent left and 
right eigenvectors r° and f are r° = c*{l + Xri) and 

= 1 + XH, and Xw, Xn and xii are uniformly small 
for a given sufficiently small e under all variations 
of rate constants. Here, w* is the stationary cycle 
reaction rate and c* are stationary concentrations 
for a cycle (fTTj) normalized by condition — 1- 

The estimate ew* for — Aq is e-small with respect to 
any reaction of the cycle: w* = kic* < Ki for all i 
(because c* < 1), and ew* <C Ki for all i. 

The kinetic equation for the perturbed system is: 

ci = -(1 + e)KlCl + K„C„, 
Cj = -KiCj + Ki_ic,j_i (for i ^ 1). 



1 1 



(96) 



In the matrix form we can write 



:^ Kc=iKo- efcie^e^^)c, 



(97) 



where Kq is the kinetic matrix for non-perturbed 
cycle. To estimate the right perturbed eigenvector 
and eigenvalue Ao we are looking for transformation 
of matrix K into the form K = — 9re^^ , where 
if is a kinetic matrix for extended reaction system 
with components Ai, ...An, K^r = and = 1. 

In that case, r is the eigenvector, and A = —6ri is 
the correspondent eigenvalue. 

To find vector r, we add to the cycle new reactions 
Ai — > Ai with rate constants eKir^ and subtract the 
correspondent kinetic terms from the perturbation 
term ee^e^^ c. After that, we get K = — Ore^^ 
with 9 = eki and 

{Krc)i = -fcici - efci(l - ri)ci + fc„c„, 
{Krc)i = —kiCi + ekiTiCi + fci_iCi_i for i > 1 

(98) 

We have to find a positive normalized solution > 
0, Xj = 1 to equation KrV = 0. This is the fixed 
point equation: for every positive normalized r there 
exists unique positive normalized steady state c* (r) : 
KrC* (r) = 0, c* > 0, X, c* (r) = 1. We have to solve 
the equation r = c* (r) . The solution exists because 
the Brauer fixed point theorem. 

If r = c*(r) then fc^r^ — efcir^ri = fci_iri_i. We 
use notation w* (r) for the correspondent station- 
ary reaction rate along the "non-perturbed route" : 
w*{r) — kiTi. In this notation, w*{r) — €riw\{r) — 
w*_-^{r). Hence, \w* {r) — wl{r)\ < £wl{r) (or Ihri — 
fciril < efciri). Assume e < 1/4 (to provide 1 — 2e < 
1/(1 ±e) < l + 2e). Finally, 



(1 + xOc* 



(99) 



where the relative errors \xi\ < 3e and c* = 
c*(0) is the normalized steady state for the non- 
perturbed system. For cycles with limitation, w 
(1 -I- Xi)kum/ki with \xi\ < 3e. For the eigenvalue 
we obtain 



Ao = -ewlir) = -ew*(r)(l + 

= ~ew*{l + x)^-ek,c*{0){l+x) 



(100) 



for all i, with \q\ < e and \x\ < 3e. |x| < 3e. 
Therefore, Aq is e-small rate constant ki of the non- 
perturbed cycle. This implies that Aq is e-small with 
respect to the real part of every non-zero eigenvalue 
of the non-perturbed kinetic matrix Kq (for given 
number of components n). For the cycles from mul- 
tiscale ensembles these eigenvalues are typically real 
and close to —ki for non-limiting rate constants, 
hence we proved for Ao even more than we need. 

Let us estimate the correspondent left eigenvector 
l'^ (a vector row). The eigenvalue is known, hence 
it is easy to do just by solution of linear equations. 
This system of n — 1 equations is: 



- li{l + e)ki+l2ki = Xoh 
~ hh + k+ih = ^ak, i = 2, ...n - 1. 
For normalization, we take = 1 and find: 



(101) 



h, k+i — { 



1 1 i>2. 
(102) 

Formulas dM]), pOO]) and (fT02)) give the backgrounds 
for surgery of cycles with outgoing reactions. The 
left eigenvector gives the slow variable: if there are 
some incomes to the cycle, then 



Ci = -(1 + e)KiCl + K„C„ + 01 (t), 

Ci = -KtCi + Ki-id-i + (j)i{t) (for i ^ 1) 
and for slow variable Z— ^liCi we get 



Ac ^ 



' k4>i{t). 



(103) 



(104) 



This is the kinetic equation for a glued cycle. In 
the leading term, all the outgoing reactions Ai ^ 
with rate constants k = eki give the same eigenvalue 
-ew* (fTOO]) . 

Of course, similar results for perturbations of zero 
eigenvalue are valid for more general ergodic chemi- 
cal reaction network with positive steady state, and 
not only for simple cycles, but for cycles we get sim- 
ple explicit estimates, and this is enough for our 
goals. 
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